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ABSTRACT 

The GEOSTAR-I multiple arc geopotential and station 
position estimation system is described. A detailed pres- 
entation of the mathematical model and formal documenta- 
tion of the principal program components are included. 
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GEOSTAR- 1 

A GEOPOTENTIAL AND STATION POSITION RECOVERY SYSTEM 


by 

C. E. Velez and G. P. Brodsky 
Goddard Space Flight Center 


I. INTRODUCTION 

The GEOSTAR- 1 system is a multiple arc, multiple satellite geopotential coefficient and station 
position recovery system. Its principal feature is the capability to process the full spectrum of 
available tracking data for the determination of a number of geodetic parameters, simultaneously 
utilizing the accuracy of the measurements and maintaining computational efficiency. 

The GEOSTAR-I system is basically comprised of: 

(1) A modified version of NONAME which is a single arc, definitive orbit and station recovery 
system consisting of ODP and ancillary data handling programs (Reference 1). 

(2) The matrix algebra programs MERGE, SOLVE and EIGENVALUE from the multiple arc 
lunar potential estimation program LUNGFISH (Reference 2). 

GEOSTAR-I is designed to operate in either a single or multiple arc mode, hi the single arc 
mode, the parameter set is solved for by using a differential correction process with a least squares 
estimator which uses the a priori estimates of the parameters and their covariance matrices. 

In the multiple arc mode, the single arc program generates a system of normal equations which 
corresponds to each individual arc and parameter set. The matrix algebra programs combine and 
process these systems of equations for a multiple arc, linear least squares solution, or essentially 
one differential correction iteration. 

Under the current capabilities of the system, each arc can have a parameter set consisting of 
a maximum of 50 parameters of the following types: 

• Geopotential coefficients 

• Tracking station coordinates 

• Arc dependent parameters including state, drag and solar constants, and, in the single arc 
mode, tracking instrument errors. 



The maximum number of parameters which can be processed in a multiple arc solution is 500. 
Moreover, by using the same system of normal equations, the determination of any subset of the 
original parameter set is permitted by the "suppress" feature of the system. 

Several optional methods, available in the system, enable the detection and analysis of an ill- 
conditioned normal matrix. These methods include matrix pseudoinversion with or without rank 
reduction, gradient, or steepest descent methods, and parameter transformations in which, by a 
canonical decomposition of the normal equations, the linear combination of parameters well- 
determined by the data can be detected. 

Some of the principal capabilities of the single arc program include: 

• Cowell type numerical integration of the equations of motion and linear variational equations 
in which the integrator used for position and velocity is independent of the integrator used 

to obtain the partials of position with respect to the parameters. 

• Simulated data (without noise) generation, and processing capabilities, including the proc- 
essing of rectangular coordinate data. 

• Complete flexibility in the choice of parameters to be estimated, so that in the single or 
multiple arc operations, arc dependent parameters may be excluded. 

• Position partials generation capability in the orbit generator mode of the program, permit- 
ting independent investigations on the sensitivity of the satellite position with respect to 
variations in the geopotential harmonics. 

• All the current data reduction and analysis capabilities of the single arc operational NONAME 
system. 

• An optional variable stepsize integrator, providing efficient integration of the equations of 
motion and the variational equations associated with satellites having high eccentricities. 

The GEOSTAR-I system, written in 360 FORTRAN, is currently operational on the IBM 360, 
Models 95 and 91. The operating instructions, system testing procedures and test results can be 
found in Reference 3. 
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II. SYSTEM DESCRIPTION 


The primary module of the GEOSTAR-I system is a modified version of the single arc orbit 
determination and geodetic parameter estimation program from the NONAME system. All the cur- 
rent data reduction capabilities, operating instructions, and overall design ox this program have 
been maintained. The principal modifications made to this program, as discussed in the following 
sections of this report, are the incorporation of improved numerical integration techniques, the com- 
putation of partial derivatives for geopotential coefficient estimation, and interfaces with LUNGFISH 
multiple arc matrix algebra programs. The details of the present single arc NONAME ODP not 
directly related to these modifications will not be discussed in this report but may be found in 
Reference 1. 

In the single arc mode, the GEOSTAR-I ODP can be used to estimate orbital and geodetic pa- 
rameters from satellite tracking data by differential corrections using a Bayesian least squares 
method of estimation. The maximum single arc parameter set permitted is 50 parameters, op- 
tionally selected from the following: 

• Geopoential coefficients through C 30 30 and S 30 30 (maximum 50) 

• Tracking station coordinates (maximum 30} 

• State, or epoch position and velocity vectors 

• Physical constants related to atmospheric drag and solar radiation forces (maximum 2) 

• Tracking errors, including measurement biases and station timing biases (maximum 44). 

In the multiple arc mode, the single arc ODP is used to construct a system of normal equa- 
tions for a specified parameter set from the given parameter estimates and observations for each 
arc. These normal equations are in turn converted to the LUNGFISH matrix processors' required 
format (B matrix), and written on tape for use in either the MERGE or SOLVE programs. Each B 
matrix tape contains arc identification, normal equations, parameter estimates, and parameter 
identification labels defining the parameter set associated with the corresponding normal equations. 

The MERGE program is used to copy the various single arc B matrix tapes ontc one logical 
merged tape which may consist of several physical reels of tape. The resulting "merged" B ma- 
trix can then be input to the SOLVE program for a multiple arc solution. 

More precisely, the SOLVE program performs the following functions: 

(i) On option, selected parameters are "suppressed" or deleted from each individual in- 
put B matrix, thus providing selectivity in the parameters to be used in a particular mul- 
tiple arc solution; 

(ii) The arc parameters, if present, are then eliminated from each B matrix, producing 1) a 
set of arc independent or "reduced" matrices which are used to solve for the arc in depend- 
ent parameters; and 2) a set of "backsubstitution matrices” which are used to solve for 
the arc dependent parameters; 
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(iii) The reduced matrices are then combined (essentially added), producing the "combined" 
matrix which is then inverted to produce a solution set for the gravity and station 
parameters; 

(iv) Finally, the arc dependent parameters are computed for each arc using the saved "back- 
substitution” matrices and the arc independent solution set. 

The SOLVE program essentially completes a single iteration of a multiple arc differential cor- 
rection process. This program outputs the updated estimates of the gravity station position and 
arc parameters in both printed and punched card form. The card output is in the GEOSTAR-I ODP 
format in order that successive iterations may be readily performed. Optional output includes the 
variance-covariance and correlation matrices for the parameters as well as any of the interim 
matrices used in the solution, including tapes containing the combined and backsubstitution matrices. 
These tapes can then be used directly in SOLVE with an additional B matrix tape to obtain a new 
arc dependent and independent parameter set solution reflecting the arcs present on both the com- 
bined and B matrix tapes. 

An example of the above operations is displayed in Figure 1, in which arcs 1 and 2 are proc- 
essed through the MERGE-SOLVE loop, followed by a SOLVE operation using combined and back- 
substitution matrices together with an additional arc 3. 

Because of the large numbers of geodetic parameters which can be considered in a given mul- 
tiple arc solution, it may frequently happen that the resulting combined matrix is ill-conditioned 
due to high correlations or poor observability of the particular parameter sets selected. 

Various investigators have suggested different approaches to either avoid or cope with this 
problem. One of these is to determine, on the basis of an analytic perturbation analysis on the 
mean orbital elements associated with a particular arc, the magnitude of the perturbation resulting 
from a particular harmonic. A GSFC program called HAP (Reference 4) is currently available for 
this purpose. This method, however, does not currently determine the cumulative effect of linear 
combinations of harmonics, which may lead to problems resulting from unmodeled parameters. 
Other methods involving either the examination of the correlation matrices or the use of con- 
straining a priori information can also be useful tools for this problem. 

In addition to these methods, algorithms have been developed (References 5 and 6) which, in 
certain cases, enable one to either precisely determine which parameters can be estimated with 
physical significance by direct examination of the normal system, or avoid numerical difficulties 
by obtaining least square solutions which do not involve the inversion of a poorly conditioned 
matrix. Several well-known methods of this type have been made available in the GEOSTAR-I 
system. 

The first method involves the diagonalization of the normal matrices and uses the eigenvalues 
and eigenvectors to determine the parameters or linear combinations of parameters which can be 
estimated with physical significance. The canonical decomposition of the normal matrix into eigen- 
vectors and eigenvalues is performed by the program EIGENVALUE which accepts the combined 
matrix from SOLVE as input. The EIGENVALUE program will eliminate any station or arc 
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The second method involves the computation of a unique minimal norm least squares solution 
of the normal equations by using a generalized inverse in lieu of a Gauss-Jordan inverse of an ill- 
conditioned matrix. The generalized inverse used is the Penrose pseudoinverse and is computed 
using the Andree Algorithm. Details concerning this method can be found in Reference 6. This 
option has been incorporated into the SOLVE program so that one can elect to use either the Gauss- 
Jordan elimination subroutine for a direct inverse, or the Andree subroutine for the pseudoinverse 
of the combined matrix. 
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Finally, in the single arc ODP, an option is available to use a gradient or "steepest descent" 
method whenever an ill-conditioned matrix is encountered during the iterative process. This 
method is given in detail in Reference 7 and outlined in Section 3.1 of this report. 

Before proceeding with the description of the mathematical methods used in the GEOSTAR-I 
system, it is noted that future GEOSTAR systems which are currently under development or plan- 
ning will include the following capabilities : 


• A maximum of 250 parameters for a single arc 

• The use of a priori information in the multiple arc differential correction process 

• The incorporation of analytic partial derivatives of the elements with respect to the geopc- 
tentiai coefficients 

• Techniques to handle pass dependent instrumentation biases efficiently in a multiple arc 
environment 

• The use of advanced numerical integration techniques, such as generalized multistep 
methods and multistep starters, to improve the efficiency of the numerical integration 
process 

• The use of semi-analytic techniques to solve the linear variational equations in which the 
state transition matrix is utilized in the solution for the partials. 

• The use of advanced data management methods to handle the large amounts of tracking data 
required for a full geopotential solution. 
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m. METHOD 


3.1 Geopotential Coefficient Estimation 

Essentially, the geopotential coefficients are estimated by augmenting the existing NONAME 
ODP normal equations by including these parameters. This augmented system of normal equations 
is solved within either the framework of the current single arc differential correction process in 
NONAME, or the multiple arc processing within the LUNGFISH system. 

The normal equations in the GEOSTAR-I ODP are of the form: 

[a t WA + P 0 _1 ] Ax = A T WAO - P 0 _1 ( 2 Ax) ( 1 ) 

is the number of observations, and n the number of parameters, then 

an m xn matrix of measurement partials with respect to the parameters to be estimated, 
i.e,, A = {a. . } , where, if m. is the i th observation and x. is the j th ^ parameter to be 
estimated, then 

i = 1 , 2, • • • m 

a ij ~ <9x. 

1 3 = 1 , 2 , » 

W = matrix of measurement weights 

P 0 = a priori covariance matrix on the a priori estimates of the parameters 
Ax = n x 1 parameter correction vector to be determined 
Ao = m x 1 residual vector 

2 Ax = the accumulated correction vectors over previous iterations. 

These normal equations are precisely those which exist in the NONAME ODP with the excep- 
tion of the quantities needed to estimate the geopotential coefficient parameters C n m and S n m . 
Assuming that the a priori estimates and covariances of these parameters are available, the addi- 
tional quantities required are of the form 

3M(t) 

dx. 

1 

where M( t) is some measurement at observation time t, and x. is one of the harmonics. 

By the chain rule, 


where if m 
A = 


3M(f) 

dx. 

J 


<9M(t) <?x(t) 
dx dx. 


dM(t) dx(t) 

di dx i 


( 2 ) 
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where x = (x, y, z) and x = (x, y, z) are the satellite position and velocity vectors at observation 
time t . Note that here and throughout this report, the following notation is used for derivatives 
with respect to vectors: 

If T and g are vectors, 7 = (f x , f 2 , f 3 ), g = (g t , g 2 , g 3 ) and h is a scalar function, then 


<9h _ (dh dh 5h\ 

<?¥ ” ’ dg 2 • dg 3 ) 


dg _ ( 3g l dg 2 , 

3h “ coi \ 3h - 3h ' dhj 


and 



££ i 

£1 

Hi 


5 <n 

dg 2 

ag 3 

fill . 

£1 

£1 

££a. 

3gJ " 

a *i 

ag 2 



£13 

af 3 

di 3 


dg. 

ag 2 

d % 3 


Now the instantaneous observation partials with respect to satellite position and velocity 


aM(t) 

dx 


<9M(t) 

dx 


are computed precisely as in the current NONAME ODP for all the observation types considered 
in the system. Hence by Equation 2, all that remains to be computed to form the full observation 
partials are the position and velocity partials 

ax(t) ax(t) 


The equations used to compute these position partials can be easily derived by examination of 
the equations for the two-body and perturbative accelerations acting on the satellite. These ac- 
celerations can be e~ pressed as 


au(x) 

dx 


+ F dras (X, 


. W 


A sun 


(3) 
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where U(x) is the geopotential given by 


U(x) 


30 " 
™EE 

n = 0 m = 0 



[ C n* 


cosm\ + S sin m M P m (sin 


(4) 


where 


r - radius from center of earth to satellite ' 
a c - earth's semi-major axis 
GM - gravitational constant times mass of earth 
A. - geocentric longitude (positive east) 

>// - geocentric latitude 
P“ (simp) - associated Legendre polynomials 

and where F drag , F sr , F sun , F moon are the acceleration vectors due to atmospheric drag, solar radia- 
tion pressure, solar and lunar gravity respectively. The equations for these forces are precisely 
those used in NONAME and can be found in Reference 1. 

Now if x. represents any physical parameter occurring in the right hand side of Equation 3, or 
a state vector element, then the variation of the acceleration with respect to this parameter can be 
expressed as 


dx _ 3x. dx dx dx dx 

dx. ~ dx dx. + 3-6- dw. + dx: 

i i ax i 1 

obtained by differentiation, noting that by Equation 3, x is a function of position and velocity. By 
interc han ging the order of differentiation, this equation can be rewritten as 

d 2 _ dx(t)l _ dx dx dx 

dt 2 [TvJ " d * 3x i + <?* 

which is a 2nd order linear differential equation in the position partial ax(t)/3x j i.e., if vector 
functions y (i) , y (i) , 7 ( ' 5 are defined by 


_d_ 

dt 


<9x 

dx. 


(5) 


y (i ) (t) 


dx(t) 

dx . 
i 


y (i) (t) 


dx(t) 

dx. 

j 
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f (i) (t) 


dx 

dx. 

j 


and the matrix functions A( t ), B( t ) by 


A(t) = 


B(t) = 


then by Equation 5, y (i) (t) satisfies the equation 


C?x( t ) 
dx 


d*(t) , 

dx 


5T (i, (t) = A(t)yO') (t) + B(t)yd) (t) + (t) , (6) 

which is called the linear variational equation with respect to the parameter x. . Note that the 
Jacobian matrices A(t) and B(t) are independent of the particular parameter under consideration, 
so that they are the same for all the parameters. These matrices, as currently computed in the 
NONAME program, are obtained from the differentiation of Equation 3 with respect to position and 
velocity, where only the more significant contributions are retained to improve efficiency. For 
example, the term 

dx j_dxj 


only includes tesserals to third order and zonals to fourth order; see Reference 1 for details con- 
cerning these approximations. 


The terms in Equation 6 which are not currently defined in the NONAME ODP are the forcing 
terms T (i) (t) = dx/dx. , where x. is a geopotential coefficient c n>m ors n m . These are computed 
as follows: 

Letting 4> denote the vector of satellite spherical coordinates, 4> ~ (r, \p), from Equation 3 it 

follows that 


dx 

e 

<9U 

! d 

'dV d$ 

dC 

nm 

5 c nm 

dx 

' " dC 

nm 

.a? dx _ 
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likewise, 



where d\j/d<p is the vector 


5U 

3? 


/ 3U dV 3u\ 
\ 3r ’ d\ ’ d4>) 


and (3^/ax) is the transformation matrix given by 


d r 

dy 

dk 

dy 

d\)j 

dy 


where 


s = (x^y 2 ) 1 ' 2 , 


and where, by differentiating Equation 4, the partials 


d 


i d 

dU 

d 

3U 

d 

3U 

dC 

nm 

\d<p/ 

V^c nm 

dr ’ 

dC 

nm 

d\. ’ 

dC 

nm 

d\p 

d 

(dU\ 

( d 

3U 

d 

3U 

d 

3U 

dS 

nm 

£ 

i 


dr 5 

d 2 

n m 

dk ' 

3s nm 

nm 

d\ p 
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are given by: 


3 

3U 

Jc 

dv ~ 

nm 

d 

(9U 

t?C dK 

nm 

d 



3C W 

nm r 


3 

5U 

3S 

n m 

<9r 

3 

3U 

3S nm 

n m 

dk ~ 

8 

au 

3S nm 

n m 

Si/' 



(n + 1) cos mA. P n m (sim/ 1 ) 



n 

(m sin mA.)P n m (sin 


cos mA. 



(sin i p) - m tan i^P n m (sin i/») 





n 

(n + l)sinm\ P n m (sin \ ft) 


(m cos mA.) P^ 1 ( sin i/') 


sin mA. 



(sini p) - mtani/iP" (sini|i) 


(10) 


and the forcing functions are then computed using Equations 7, 8 and 9. 

Using these forcing terms in Equation 6, the required partials with respect to geopotential 
coefficients as well as other parameters are obtained by numerically integrating this equation for 
each j , using the initial conditions 


y (i) (t 0 ) = ^ (i) (t 0 ) = 0 

for all j corresponding to geopotential coefficients. 

The measurement partials are then formed using Equation 2 and used in the normal system 
Equation 1. Upon completing the formation of these equations, the ODP will either output them on 
a tape in the SOLVE format for a multiple arc solution, or perform an iterative single arc solution 
for the parameter correction vector Ax. In the single arc process, the normal matrix 

S = ^ A T WA + P Q - 1 J 


is first inverted using the Gauss-Jordan method with partial pivoting. The resulting matrix is 
tested to determine the condition of the matrix S. If it is found to be ill-conditioned, the solution 
is determined by the gradient technique. This is essentially a step in the method of steepest 
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descent for obtaining least square solutions , which finds the scale-factor multiple of the right hand 
side of Equation 1, denoted by b, which minimizes the quadratic form 

R(Ax) - (SAx-b) T (SAx-b) 

i.e., a scalar X is found such that if Ax = Xb , then R(Xb) is minimized, or 


This solution is found to be 


x 


b T ■ 


T 

y y 


where 


y - Sb , 

and the correction vector is then given by Ax = Xb . 


In either case, the parameter set is then updated and the process is repeated until a convergence 
criterion is met, or a maximum iteration count is exceeded. 

The multiple arc process is described in Section 3.3. 

3.2 Numerical Integration and Interpolation 

The numerical integration of the equations of motion (3) for position and the variational equa- 
tions (6) for the position partials, is performed by using a summed Cowell type integration method. 
These formulas are used in Lagrangian or ordinate form, and the coefficients are available in the 
system for orders 4 through 15. The starting procedure used is an 8th order Runge-Kutta, and the 
Hermite interpolation formula, utilizing functional derivatives, is used to interpolate position and 
velocity vectors as required by the process. 

3.2.1 The Integration of the Equations of Motion 

The basic integration formulas for the integration of Equation 3 are given by 



(Stormer Predictor) 


(11a) 
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X n + 1 = ll e 


II§ n + 2 ^ + 


i “0 


(Cowell Corrector) 


(lib) 


X n + 1 = h o 


o 

r s + r ° x 

" Z_< 1 


(Adams Predictor) 


(lie) 


x. , = h 

n + 1 e 


3 + y^ft 

3 n / | i n + 1 - i 


(Moulton Corrector) 


(lid) 


where 

h D = the stepsize used for the integration of the equations of motion, 
k e = p - 2, where p is the order of the formulas (11), 
r S , '“S = the first and second sums of the accelerations, which can be defined as 

n * n 7 

X S = V " 1 5 ? 

n n 

n S = V ' 1 J S = V“ 2 x ’ 

n n n 

a. , a.*, /S. , j3* = the ordinate integration coefficients which depend on the order p. 

The ordinate form coefficients can be formulated by defining the difference form coefficients 
recursively and converting to ordinate form. In the GEOSTAR-I system, these coefficients have 
been precomputed in rational form (Reference 8) and entered into the system as a permanent data 
file for orders 4 through 15. 

The algorithm used to integrate the equations of motion is the following: 

(i) Compute a set of "starting" values for the accelerations 

x k -i • 1 = l. 2 . •“ k e • 

e 

and sums l s k , n s k using an independent procedure; 

(ii) Using Equations 11a and 11c with n = k e , obtain predicted values 


X W 
n+ 1 


;( p ) ; 
n+ 1 ’ 


(iii) Evaluate the force model using the last computed position and velocity vector to obtain 


n+1 5 
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(iv) Using Equations lib and lid, obtain corrected values 


(v) Compare the magnitude of the vector 


X 

n + 1 


•<p> 


with a predictor-corrector tolerance. If this difference vector is sufficiently small, the 
predictor -corrector cycle is complete and step (vi) is then executed. If it is not suf- 
ficiently small, the predicted values are replaced by the corrected values and steps (iii) 
to (v) are repeated. The maximum number of iterations allowed is three. Note that it is 
possible to complete a predictor- corrector cycle with only one force model evaluation. 

(vi) Compute updated sums by 


IS n + 1 

- X , . 
n + 1 

+ 

n s n+1 

= IF 

n + 1 

+ n s 

n 


completing the integration step. Steps (ii) through (vi) are then repeated with 

n = n + 1. 

3.2.2 Integration of the Linear Variational Equations 

The integration formulas used to integrate the linear variational Equations 6 are given by 


V (i) 

y n + l 


yj-ll ~ h v 


v 

npO') + 2]a i *y n ( i)- i 

i “0 
k 

v 


(12a) 


j = 1, 2, 


(12b) 


L i-o J 

where indicates the j th position partial at t n+l , and 

h v = the stepsize used for the integration of the variational Equations 6, 

k v = p - 2, where p is the order of the formulas (12), 

I P n (j ) . n P n (i) = the first and second sums of the acceleration partials with respect to the j th 
parameter, 
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a.*, >i* = the ordinate corrector integration coefficients as in formulas (lib) and (lid), de- 
pending on the order of formulas (12). 

The GEOSTAR- 1 system was designed so that the stepsize and order used to integrate the vari- 
ational equations need not be the same as those used to integrate the variational equations. This 
allows efficient utilization of the varying accuracy requirements between the position and position 
partials. For example, in certain cases, it may be possible to obtain sufficiently accurate position 
partials using two or three times the stepsize used for the satellite position, or half the order, 
thereby improving overall efficiency. 

The method used to integrate the variational equations employs a closed form solution of Equa- 
tions 12, or the "corrector-only" technique. The required equations are derived as follows: 

From Equation 6, a typical variational equation at time t n+1 can be expressed as 


y n + l = A n + 1 y n + l + B n + 1 y n + l + f n + l' 


(13) 


where for the moment, the j superscript is dropped. From Equations 12 it follows that 


Letting 


y . h 2 

J n + 1 v 


y ,, - h 

J n+1 v 


v 

ir P + a * v + / a * v 

n 0 y n + 1 / , i y n + l-i 


v 

r P + & * y + T B * V 

r n ^ n + 1 / y n + l-i 


X = h 2 


V„ = h 


up *** / ' a..* v + cl * f 

n / , i y n+l-i 0 n+1 

i=l 


X P + / B.* v . + (3 * f 

r n / t y n + l-i H 0 n+1 
i“l 

and substituting Equation 13 into the right hand side, it is possible to rewrite Equation 14 as 


( 14 ) 


(15) 


>n + l 


h v 2a 0* A n + l y n + l + H v 2 < B n + I y nn + X n 


( 16 ) 


'n + 1 


h v K A n+1 y n + 1 + h v AT B n + 1 y n + l + V n 
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which can be expressed as the matrix equation 


y'ntl 

_ 

h v 2a o*A.n 

h , 2 *0*B n + 1 



+ 

X 

n 

^n + 1 


f 'o A r. + 1 

h v ‘*0 ^n + 1 


y„ + i 


V 

n 

— — 


— 







or letting h be the 6 x 6 matrix 


H 


\WB n+1 

h v @0 ^n + 1 h v ^0 B n + 1 


and I the identity 6x6 matrix, then the 6 x 1 vector 



is the solution of the equation 


and hence is given by 


_ — 




y n+ i 


X 

n 

• 



1 


V 

n 




. — — 



- [I _ H] -1 

X 

n 





( 17 ) 


(18) 


(19) 


We note that the matrix H is independent of the particular parameter being considered so that 
the inversion in Equation 19 is performed only once per step regardless of the number of parameters 
being estimated. Using Equation 19, the solution of Equations 12 for all the position, and velocity 
partials can be expressed as 


y CD 
y n + l 

v <?> • • 

y n + l 

• 


XCD 

n 

X CD • • 
n 

..XCN> 

V 

’ y n + l 

t (2) .. 

y n + l 


6XN 

V CD 

n 

V C2) . . 
n 

• • V CN) 

n 


( 20 ) 


6XN 
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where N is the total number of partials being computed, and where the x n (i) and v n (i> are defined in 
the natural way from Equations 15 as 


xu> 

n 


K“> 

i ~i 


Wl - i 


+ a * f 


' <i) 

n + 1 


(21) 


V <i> 
n 


h 


V 





(i) 

n + lj 


The algorithm used to integrate the system of n linear variational equations is the following: 
(i) Compute a set of "starting" values for the acceleration partials 



i = l, 2, • • • k 

V 


and sums n P k (i) , using an independent procedure, for each j = 1, 2, • • • N. 

(ii) Letting n = k compute the values of satellite position and velocity at time t = t. +h , 

K V V 

denoted by x n ' +1 and x n ' +1 . Normally, if h v = h e , these values can be obtained by integrating 
the equations of motion as described in Section 3.2.1 to time t = t n+1 , i.e., x n ' +1 = x n+1 and 
*„Vi “ ¥ n+1 . If h„ / h e , these values are obtained by interpolation, where the integration 
of the equations of motion continues until enough position and velocity data is available. 

(ill) Using the vectors x n ' n , x n ' +1 compute the Jacobian matrices A n+1J B n + 1 and the forcing terms 

7 (i) • 

r, + l > 

(iv) Using the mati’ices A n+1 , B n+1 , form [i -h] by Equation 17 and using the vectors 


f W 

n+1 


y n ( j ) 


ip (i) 


np (j) 

n 


compute the vector quantities 

X <i) , V<i> j = 1, 2, ••• N 

n n * 

using Equations 21. 

(v) Invert [l -h] and perform the multiplication in Equation 20 to obtain the values for the N 
position and velocity partials, completing the integration step to time t n+1 . Steps (ii) 
through (v) are then repeated with n = n +1. 

3.2.3 The Integration for the Case of Velocity- Free Accelerations 

For the case in which the accelerations acting on the satellite do not depend on the velocity 
vector, i.e., F drag in Equation 3 is zero, the number of computations performed for the integration, 
as described in Sections 3.2,1 and 3.2.2 can be considerably reduced. 
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In this case, for the integration of the equations of motion, Equation 11c which predicts the ve- 
locity vector is not used, and Equation lid, the velocity corrector, is used only once, after conver- 
gence of the position corrector formula Equation 11c. 

For the integration of the variational equations, note that the Jacobian matrix B(t) in Equa- 
tion 6 is zero, so that Equation 13 is simply of the form 

y n + l “ A n + 1 y n + l + ^n + 1 ’ 

By using the notation of Section 3.2.2, the corrector formula to be solved can be expressed as 

y n + l = h v 2 a o’ A „ + l y n+l + X n 

with solution 

y . . = [I -H]' 1 X - 

where H is now the 3x3 matrix h v 2 a Q * A n+l , and I the 3x3 identity. The velocity components are 
obtained from 

y n + i = K AT A n + l y n + l + ’ 


Hence, for this case, Equation 20 is replaced by the equations 


3xn = 

7(f) = fh /3 * A 1 [V <*> ••• v < N >1 + [><« V< : 

_ y „ 1 y " 1 yn + 1 J 3 XN [ v 0 n + 1 J 3 X3 [ n+1 yn + 1 J 3 XN L " " 


which can be seen to require significantly fewer computations. 


3.2.4 The Starting Pro cedur e 

■ The method used to compute the required starting arrays of accelerations and acceleration 
partiais is an eight order, single-step Runge-Kutta type numerical integration. 

This method is defined as follows : 

Given an arbitrary system of first order differential equations 


x = f(t, x) 


(23) 
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with initial vector x(t 0 ) = x 0 , the solution vector at time t x = t 0 +h t , an arbitrary stepsize, is 
approximated by 


9 

*(* 1 ) = *(?o) + h i Ci " k " i ’ 

i^O 


where, for each i, 


(24) 



i = 1, 2, ••• 9 


(25) 


and where a., b. . and c. are constants defining the method. These can be found in Reference 9 
and are also given in the subroutine description RK in Section 4.2. 

Note that to compute x(t 2 ) = x(t x +h 2 ) where h x r h 2 , it is possible to repeat the above process 
starting with x(t t ) and using h 2 in the definition of the k. (Equation 25, so that the method ad- 
vances from step to step independently of the stepsizes taken, as is usually the case for single- 
step methods. 

Now, by the usual procedure, the three-dimensional system of 2nd order equations given by 
Equation 3 can be expressed as a six-dimensional system of 1st order equations; likewise, the var- 
iational Equations 6 can be reduced to a system of N 1st order equations of dimension 6, so that both 
systems (Equations 6 and 3) can be reduced to a system of the form given by Equation 23. 

The algorithm used to form the starting values for the multistep integrations is the 
following: 


Case I: h = h 

v e 

In this case, the number of timp points produced by the single step method is l - max (k v , k e ) , 
i.e., using the given initial state for the equations of motion and the initial values 

y <j) (t 0 ) - = 0 . 

for the variational equations, the single step integrator is used to compute the values: 


X 



-at t. 


y .(», 


y U) , y.U) 


t. + ih 

0 e 


i = 1 , 2 , ••• l 

j = 1, 2, • • • N * 
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Case II: h v f h e 

In tMs case, the number of time points produced by the single step method is l - k v + k e , i.e., 
the single step integrator is used to compute the values: 

X. , X. , X. at at t. = t 0 + ih e , i, = 1, 2, • • • k e 

y. (i) , y.(i>, y.<i) at t. = t„ + ih , i = 1, 2, k 

j = 1, 2, • • • N . 


In either case, the stepsize used by the single step method is independent of the stepsizes h v or h e , 
and is generally chosen to be a fraction of the stepsize, h e ,used to integrate the equations of motion. 

Once the above starting values have been obtained, the starting first and second sums are 
computed as follows: 

First, the sums 


Ic Ip (i ) lip 

°k -1’ r k -1’ I 


(j) 


are computed by inverting the respective corrector formulas (lid), (lib), (12b), and (12a), with 
n = k e - 1 or n = k y - 1 . 

For example, from Equation lid, 



i=0 


where we see that all the quantities on the right side are given by the starting values. Once these 
have been computed, the required starting sums are given by 




I S„ 


J k -1 


and likewise for I P 1 < 3> , n P v (i) . These equations complete the starting procedure. 

*v v 


(26) 


3.2.5 Variable Stepsize Integration 

A version of the GEOSTAR-I ODP program is available which uses an integration method 
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allowing dynamic stepsize modification during the integration. This version is designed specifically 
for the integration of the equations of motion and the concomitant variational equations associated 
with high eccentricity satellites, where normally fixed stepsize integration methods are extremely 
inefficient. 

In the variable step integration program, the equations of motion and variational equations are 
integrated with the same integration methods outlined in Sections 3.1.1 and 3.1.2, with the exception 
that h e = h v = h , i.e., a common stepsize is used. 

The stepsize variations axe based on the concept of local error control, where the stepsize is 
selected so that the local error, denoted by Rn, at each step of the integration satisfies the con- 
straint equation 

t 2 < Rn < Tj (27) 

where T t , t 2 are specified upper and lower bounds. 

In the GEOSTAR-I ODP, the local error Rn is estimated by the quantity 

Rn = Cjx - x c j 

where c is an error constant depending on the order of the formulas (Equations 11) being used, and 
x n (P) , x n cc> are the predicted and finally accepted satellite position vectors respectively, computed at 
time t . 

n 

The variable step integration algorithm is then the following: 

At each step n } the test (Equation 27) is performed. We have three cases: 

(i) Rn >t x ; the stepsize is decreased, the n th computed point is rejected and recomputed with 
the new stepsize, where the required back values in Equations 11 or 12 are obtained by 
interpolation. 

(ii) Rn <T 2 ; the stepsize is increased, the n th computed point is accepted and the integration 

proceeds with the new stepsize, where the required back values are obtained by interpolation. 

(iii) Rn satisfies Equation 27; the integrator proceeds to the next step using the same stepsize. 

In either case (i) or (ii), the new stepsize is computed using the formula 



where t 3 is a specified value for the "allowable" local error :T 2 < T 3 < T : . 
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The starting procedure used for this integration method is the same as that described in Sec- 
tion 3.2.4 with the exception that the initial stepsize is modified if necessary so that Equation 27 is 
satisfied at the first Cowell step. 

Because of ca^e storage limitations, the following restrictions were imposed on this version 

of the ODP: 

(i) The maximum number of position partials allowed is 20; 

(ii) As indicated above, the stepsizes used by the equations of motion and the variational equa- 
tions are equal, although different orders may be used. 

3.2.6 Interpolation 

The position and position partials are computed at non-step points as required by the observa- 
tions, or in the integration process, using a 5 point Hermite interpolation scheme. 

Each component of the interpolated vectors is obtained by the formula 


where 


y(t) 


2^ h j ( t )y i + bj (t)y. , 
i =0 


h j (t) = 1 - ( tj )^(t) 

h. (t) = (t-t.)4. 2 (t) 


(t) 


= n 

i = n 




(ft t ) 

(^ 1 ) 



and where k = 4. 

3.3 Multiple Arc Processing 

GEOSTAR-I multiple arc solutions are obtained by using the MERGE and SOLVE programs 
with essentially the same algebraic and data manipulation techniques as used in the current 
LUNGFISH system. The ODP interfaces with these programs, their basic methods and capabil- 
ities, as well as the methods of parameter transformation and pseudoinversion, are reviewed in 
the following sections. 
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3.3.1 ODP Interfaces; Matrix and Parameter Set Handling' T? ? ocedures 

The GEOSTAR-I ODP— MERGE— SOLVE interfaces can be summarized as follows: 

(i) B Matrix— The principal data link between the ODP program and the multiple arc proc- 
essors MERGE and SOLVE. A B matrix is defined as a matrix containing the known 
elements of the normal equations (Equation 1, Section 3.1), i.e., the B matrix contains the 
elements of the matrix 

s = [a t wa + p 0 -1 ] ’ 

and the right hand side vector 

b = A T WAO - P 0 _1 (2 Ax') . 

The elements of the B matrix are arranged so that the geopotential and station position- 
parameters precede the arc parameters. Figure 2 displays the parameter sequences as 
they occur in the ODP normal matrix and as they are required in the B matrix. In ad- 
dition to the normal equations, the B matrix tape contains arc and parameter identifica- 
tion labels, as well as the parameter values. This labeling scheme is also indicated in 
Figure 2. The precise format specifications for this tape are given in Appendix B, with 
further description given in the subroutine WTBMAT in Section 4.2. 

(ii) SOLVE Punched Card Output— After a multiple arc solution is obtained, the SOLVE pro- 
gram will output a deck of punched cards containing the updated arc dependent and in- 
dependent parameters. These cards are in the input format required by the ODP program 
for the parameter estimates, facilitating an iterative process. It is noted that although 
SOLVE performs the station location calculations in rectangular coordinates (Figure 2), 
the updated station location parameters are converted to spherical coordinates for output. 

The essential matrix and parameter set handling procedures available in the GEOSTAR-I 
system are: 

(i) Merging— The MERGE program is used to copy B matrices and parameter set matrices 
from up to four input tapes onto one tape. The output tape used may already contain 
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B matrices, in which case the input 
matrices are copied after them. Also, 
the input tapes may contain several 
B matrices each. The "merged" tape 
can then be input into the SOLVE pro- 
gram, where the arcs used in the solu- 
tion can be any subset of the total set 
of arcs contained on the tape. 

(ii) Combining— The SOLVE program can 
optionally perform the following ma- 
trix handling procedures: 

(a) A combined matrix, representing 
the normal equations andparam- 
eter sets for the arc independent 
parameters over all the arcs used 
in the solution can be output in 
the B matrix format. This tape 
can then be used, together with 
an additional B matrix tape con- 
taining an arbitrary number of 
additional B matrices in a sub- 
sequent SOLVE run. 
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f "xxxx," "yyyy/" "zzzz" represent 4 digit station 
identification numbers used in the ODP. 


Figure 2— GEOSTAR ODP and B Matrix 
Parameter Sequences. 


(b) A combined matrix, as in (a), parameter set matrices for the arc dependent param- 
eters, and the "backsubstitution" matrices, over all the arcs used in the solution 
can be output on three tapes. As in (a), these tapes can then be used, together with 
another B matrix tape, in a subsequent SOLVE run, allowing the arc dependent 
parameters of a previous SOLVE run to reflect the solution of the arc independent 
parameters over all the arcs used. 

(c) Matrices containing data over the same period for the same satellite can be com- 
bined in such a way as to include the arc dependent parameters. This matrix can 
then be output on tape, as well as used to form a single arc solution. This is called 
the "COMBINE ARCS" feature. 


(iii) Multiple Reel Processing— The IBM 360 system multiple reel volume capability can be 
used with the MERGE, SOLVE programs to effectively handle large numbers of large 
matrices. The multireel capability essentially allows that tapes generated as single 
reels may be treated as a multiple reel volume for subsequent program inputs. To 
process a large number of matrices, the procedure recommended is to first use the 
MERGE program to collect as many matrices as possible onto a two or three reel 
volume, then use several of these multireel volumes as one large multireel volume 
in SOLVE. 
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3.3.2 Algebraic Matrix Operations 

The matrix operations used to determine the solution of the parameters in a multiple arc en- 
vironment are: 

(i) Elimination of Arc Parameters 

(ii) Combining of Normal Equations 

(iii) Backsubstitution 

(iv) Suppression of Parameters 

To illustrate these operations, let 

S Ax = b n = 1, 2, ••• k (1) 

n n n \*/ 

be k sets of normal equations, where the parameter sets x n are comprised of arc independent (I) 
and arc dependent (D) parameters, i.e,, 

X n = (V' X n°) n = 1, 2, k • 

where the {x n *} may or may not have common elements. For each n, we can express these sys- 
tems in partitioned form: 


\ 

c n T 


Ax 1 

n 


[b r 1 

n 

c 

n 

i 

c 

a 


Ax D 

n 


b D 

n 


( 2 ) 


To solve Equation 1 for the parameter corrections Ax n simultaneously, without forming a con- 
siderably larger matrix than any of the s n , the following procedure is used: 

First, the Elimination of the Arc Parameters operation is used to form a reduced system of 
equations. From Equation 2, we have the component equations 

(3a) 

n = 1, 2, • • • k . 

(3b) 

Solving Equation 3b for Ax n D , we get 

Ax n D = d ; 1 k D -c n AxO. (4 ) 


B Ax 1 + C T Ax D = 

n n n n 

C Ax 1 + D Ax D = 

n n n n 
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Substituting Equation 4 into Equation 3 a and rearranging, we obtain the equation in Ax n * : 

Fb -C t D~ 1 c"|Ax i = b I -C T D" l b D . (5) 

Innnnjn n nnn ' / 

Defining the matrices and vectors 


s' = b - c t d _1 c 

n n n n n 


b ' = b 1 - C T D -l b D 

n n nnn 


we obtain the system of reduced normal equations 


S' Ax 1 = b ' n = 1, 2, • • • k ’ 

n n n 11 

where each equation contains only arc independent parameters, i.e., only geopotential or station 
position parameters. 

Next, the Combining of Normal Equations operation is used to "combine" or form a matrix 
and vector union of the s ' , b ' denoted by 

n 7 n * 


s' 


i=l 


(Combined Matrix) 


b' = U b.' 

i=l 1 


(Combined right hand sides) 


which can be regarded as matrix addition, where the rows and columns of s' and elements of b' 
correspond to all the distinct arc independent parameters occurring over all the S.', say x 1 , 
Hence the elements of the S.' and b.' corresponding to the same parameter, as i varies, are 
simply added to form single elements in S' and b' . Note that the dimension of S' is then equal to 
the total number of distinct arc independent parameters. We then have the "combined" normal 
equation 


S' Ax 1 = b' ■ (6) 

where Ax 1 is the correction vector, to be determined, containing the corrections to all the distinct 
geopotential and station position parameters contained in the k equations (Equation 1). This 
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correction vector is obtained by matrix inversion 


Ax 1 = S' -1 b' 


or on option, by pseudoinversion (Section 3.3.4). 

The solution for the arc dependent parameters is then obtained by the Baeksubs titution operation: 

Ax D = D _l [b D -C Ax H ■ 

n n L n n n_j 

where each x n r is a subset of the total arc independent correction vector Ax 1 . 

The Suppression of Parameters operation is used whenever it is desirable to examine the ef- 
fects of suppressing the corrections to certain parameters occurring in either the S. matrices or 
in the combined matrix s'. This is performed by simply striking out all the rows and columns of 
the s., b. or s', b' corresponding to these parameters specified for suppression. 

3.3.3 EIGENVALUE and Ill-Conditioned Systems 
When the system of normal equations 


S Ax - b 


( 1 ) 


is ill-conditioned due to poor observability of the parameter set x selected, numerical difficulties 
may prevent a meaningful solution. A method is available which can, in certain cases, allow one 
to determine a subset of the original parameter set which is well-determined by the data con- 
tributing to s. This method is described in References 5 and 10, and is now outlined. 

The basic idea involves a coordinate transformation into a coordinate system with basis ele- 
ments consisting of eigenvectors of the coefficient matrix s in Equation 1. Let {a..} ? =1 be the 
eigenvalues of s, where sisnxn, and let {e.} " =1 be the associated eigenvectors, i.e., 


S e. - e. for all i - 1, 2, • • • n . 


( 2 ) 


Since S is a positive definite or semi-definite symmetric matrix, its eigenvalues are real non- 
negative numbers, and the eigenvectors form an orthogonal basis. Let E be the matrix 

E = [ej, e 2 . •••, ej . 

Then E is an orthogonal matrix, so that E _1 = E T , or 

e t e = ee t = i. (3) 
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Moreover, from Equation 2, E can be seen to satisfy the matrix equation 


SE = DE , (4) 

where D = diag , \ 2 , • • • \ n ). We can consider D as the matrix S expressed in the coordinate sys- 
tem with basis elements {ej? =1 ; in fact, we have from Equation 4 

E T SE = E T DE = E T ED = D . (5) 

We transform coordinates of the parameter correction and right hand side vectors by 

Ay = E t Ax 
c = E t b 


i.e., Ay and are the correction vectors Ax and b in the new coordinate system. We then have 
(since E T = e -1 ) 

E Ay = Ax 

and 


Ec - b . 


Therefore, 


S Ax - SE Ay - b 


by Equation 1, and 


E T SEAy = DAy = E T b = c 

and hence in the new coordinate system, the transformed parameter correction vector Ay satisfies 
the equation 


DAy = c. (6) 

Note however, that Dis a diagonal matrix, so that the elements of the transformed parameter cor- 
rection vector Ay which are not well-determined can be discerned by inspection. For- example, if 
Ay = (Ay 1 , Ay 2 , • • • Ay n ) , and if the i th diagonal element of D, or the i th eigenvalue of s, were 
smaller than some minimum value, the i th component of the correction vector Ay , namely A yi , 
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would not be well-determined. But note that each element of Ay is simply a linear combination of 
elements of the original correction vector Ax, i.e., 

Ay; = 


so that the above scheme can be used to establish which linear combinations of original parameter 
corrections are not well-determined. This in turn can be used to discuss which particular parameter 
corrections are not well-determined. For example, let (>- 1} ■ \ fc ) (k< n) be the set of relatively 

large or "significant" eigenvalues of S, with corresponding eigenvectors (e 1} e 2 , • • • e k ). If it is 
found that a particular parameter correction Ax. has a small coefficient in each linear combination 


e . T Ax 


e . , Ax, + e . „ Ax. + • • • + e . Ax 

x 2 i i22 inn 


n 



m= 0 


or equivalently, that the j th component of each of the eigenvectors {e.} * =1 is small, then this pa- 
rameter correction cannot be determined from the given matrix S, so that one. would then suppress 
the corresponding parameter from the parameter set, and try to solve the resulting smaller sys- 
tem. Examples of this method, as well as some possible extensions, can be found in Reference 5. 

In the GEOSTAR-I system, the calculation of the eigenvalues and eigenvectors of the real 
symmetric matrix is performed using a SHARE subroutine package called HOW. The matrix is 
first reduced to tridiagonal form using Householder's method. The eigenvalues are then computed 
using a method based on Sturm sequences, and the eigenvectors are computed using Wilkinson's 
method. 

The program EIGENVALUE accepts as input the combined matrix resulting from a multiple 
arc SOLVE run. containing geopotential and station location parameters. EIGENVALUE then elim- 
inates the station position parameters, obtaining a reduced normal matrix containing only geo- 
potential coefficient parameters. Note that this reduced matrix still contains the "effects" of the 
eliminated parameters. This matrix is then nomalized and decomposed into eigenvalues and eigen- 
vectors using the HOW subroutine package. Normalization is performed using the usual harmonic 
coefficient normalization factors 


N 

nra 


(n-ra)! (2n + 1) K/fn Tin)! -1/2 


30 



where 


Cl if m ~ 0 

K =< 

[2 if m f 0 . 

The eigenvalues and eigenvectors are then printed for user analysis of the type described in 
the first part of this section. 

3.3.4 Pseudoinversion and Ill-Conditioned Systems 

A numerical approach to the problem of finding a least squares solution of an arbitrary over- 
determined linear system 


Ax = b , (1) 

where A is airnun matrix, m > n, is the use of the Penrose pseudoinverse or a generalized inverse of 
the matrix A. Methods of this type were developed to cope with the numerical difficulties arising 
either from the formation of the normal matrix A T A, or, once formed, with the ill-conditioning or 
possible singularity of A T A . 

The Penrose pseudoinverse of A is defined as a matrix A # which satisfies the postulates: 


(i) AA # A 

= A 

(ii) A # AA # 

= A* 

(iii) (AA # ) T 

= A #T A T = A A# 

(iv) (A # A) T 

= A T A #T = A # A 


It can be shown (Reference 6) that the pseudoinverse of any matrix A always exists, is unique and 
satisfies the following: 

(i) if A' 1 exists, then A # = A -1 ; 

(ii) if 


X 


0 


A # b, 


(3) 
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then x Q is a least squares solution of Equation 1, that is, the quadratic form 


R(x) = (Ax - b) T (Ax- b) 
is minimized by the vector x 0 ; 

(iii) if more than one least squares solution of Equation 1 exists, then x 0 is the smallest in 
magnitude of all such solutions, i.e., 


for all vectors x x which are least squares solutions of Equation 1. 

The methods which have been developed to compute generalized inverses are either of the type 
where the matrix A # is formed directly from the matrix A, generally requiring more core storage 
or sophisticated algorithms (Reference 11), or, using the already formed normal matrix A T A , the 
pseudoinverse is found by the equation 


A # = (A T A) #f A T . 


In the latter case, given the normal equations 

[A T A] x = A T b = c , 

the pseudo solution would be given simply by 


x - [A T A] ^ c . (4) 

In GEOSTAR- 1, the pseudoinverse solution for the combined normal system of equations is obtained 
as in Equation 4. The algorithm used to compute the pseudoinverse of the normal matrix is called 
the Andree algorithm, which is described in detail in Reference 5 and is outlined in the subroutine 
description ANDREE in Section 4.4. This subroutine is called on option by the SOLVE program in- 
stead of the Gauss- Jordan subroutine. One of the advantages of using this subroutine is its rank 
determination capability. This feature computes the computational rank of the matrix based on a 
control parameter indicating the "noise” level of the matrix. The advantage of this feature is that 
ill- conditioning can be detected by direct examination of the matrix and the rank will be reduced 
only if necessary. If the matrix is of maximum rank, then by property (3) (i) of the pseudoinverse, 
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the pseudoinverse obtained would be equal to the inverse computed by the Gauss-Jordan method, i.'e., 
in this case 


(A t A) # = (A T A) -1 . 


If it is not of maximum rank, then numerical difficulties which would result from a direct inverse 
computation are avoided. It is remarked, however, that the applicability of the pseudoinverse to 
the problems of physical parameter estimation appears to be in question. The problem seems to 
be the determination of the physical significance of the pseudo solution in practical applications 
where noise in the data and physical correlations, as well as numerical noise, contribute to the ill- 
conditioning of the normal matrix. Studies using simulated data have been encouraging in that the 
pseudoinverse methods were found to offer distinct advantages over conventional techniques (Ref- 
erence 6). Further studies are currently underway to determine the applicability of these methods 
on problems corrupted by noise. 
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IV. NEW AND UPDATED MODULES IN NONAME ODP AND LUNGFISH MATRIX PROGRAMS 


This section describes modifications and additions made to the subroutine structures of the 
NONAME ODP and LUNGFISH matrix programs for the GEOSTAR-I system. 

To obtain the GEOSTAR-I ODP, the NONAME ODP program was modified to allow for geopo- 
tential coefficient estimation and new integration methods. To implement these capabilities, the 
following existing NONAME ODP modules were modified: 


MAIN ORBl DRAG 

COEFL PREDCT SUNGRV 

EGRAV VEVAL F(FRCS) 

ESTIM INPUT- BLOCK DATA 

DNVRT1 (DNVERT) 


and the following new modules were developed: 


For 


READGP 

SOLVGP 


Geopotential JSTORGP 


Estimation 


iSUMTOB 

IWTBMAT 


r 



ORBIT 

MMATRX 

For 

CKDIFF* 

RK 

Numerical •< 

CSTEP 

SUMS 

Integration 

ephqan 

SWTEST 


HEMINT 

TABLE 


INV2,3 
L 5 

TABLE B* 
TEST* 


To achieve the multiple arc capability of the GEOSTAR-I system, the LUNGFISH MERGE 
and SOLVE programs are used. Also, for multiple arc correlation analysis, the LUNGFISH EIGEN- 
VALUE program is used. 


The MERGE program is used without any modifications. The SOLVE program was modified to 


allow for handling of a larger number of station location parameters in rectangular and geodetic 
coordinates, to extend the "combine" matrix capabilities, to include pseudoinversion, and to satisfy 
I/O interface requirements with the GEOSTAR-I ODP. To implement these capabilities, the follow- 
ing existing SOLVE modules were modified: 


MAIN 

BEDIT 

CALTYP 

INVERT 


LBLSUP 

OPARC 

OPGRAV 

OPSTAT 


SUPRSS 
UP COMB 
ELIM 


*Subtoutines for the variable stepsize version only; see Appendix A. 
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and the following new modules were developed: 


ANDREE 

PHLINN 

OUTRAD 

The LUNGFISH EIGENVALUE program, which computes the eigenvalues and eigenvectors of the 
combined normal equations resulting from a SOLVE operation, was modified to allow for the elim- 
ination of the station parameters from the matrix before decomposition, so that the geopotential 
parameters contributing to the solution can be examined separately, while including the effects of 
the eliminated parameters. This was accomplished by adding the module ELIM. 

In the following sections, those modules listed above which are either new or significantly 
modified existing subroutines will be documented in detail. The remaining modules which received 
minor modifications are only briefly outlined, with changes indicated. Further details of these, as 
well as the unmodified NONAME ODP or LUNGFISH SOLVE subroutines can be found in References 
1 and 2. 

4.1 Modifications to Existing NONAME ODP Modules 

The modifications made to the NONAME executive program MAIN, and some of the other exist- 
ing modules, are designed to: 

• Call the new subroutines. 

« Extend data tables and variables to accommodate the geopotential coefficients through order 
30 and the new geopotential partials. 

0 Extend the position partial computation algorithms to include partials with respect to geo- 
potential coefficients. 

0 Modify the computational algorithm to allow for independent integrations of the equations of 
motion and the variational equations. 

0 Extend the program options applicable to geopotential coefficients and integration features. 

A summary of these modifications follows: 

0 MAIN — the ODP control program. Modified to control the new subroutines READGP, 

SOLVGP, STORGP, WTBMAT, INV3 and to include the I/O r equired by the 
geopotential estimation process. To allow for effective compilation, this 
program was divided into four subprograms: MAIN, OPTCRD, STATRD, 
OUTPUT. 

0 ESTIM — sums all observation partials into the normal equations matrix and right 
hand sides for each observation data point forming the normal equations. 
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• ELEMK 

COEEL 1 

• ORBl J 

• INPUT 

• PREDCT 

• DRAG 

• SUNGRV 

• FRCS 

• DNVRT1 


The subroutine was modified at entry three to exit before the solution is 
formed (subroutine SOLVGP replaces that portion of the routine), and also, 
to provide weights using a full variance-covariance matrix (using the 
VARCOV option card). 

— the subroutine ELEM which converts inertial vectors to Keplerian elements 
was modified and called ELEMK to provide the I/O in the calling sequence 
rather than COMMON. 


— modified to extend the geopotential coefficient data tables to order 30. 


— modified to extend the geopotential coefficient data tables to order 30, ex- 
tend integration coefficients (orders 4-15), extend the available program 
options, and extend data base constants requiredfor the new program options. 

— computes the observation partials using position partials. The subroutine 
was modified to accommodate the additional measurement partials with 
respect to the geopotential coefficients, and six more observation types (rec- 
tangular coordinate measurement types-x, y, z , x, y , z). 

— computes the acceleration of a satellite due to drag forces. Modified to 
not compute density above 1000 kilometers. 

— computes the acceleration on a satellite due to the gravitational attraction 
of a disturbing body. Modified to save variables for the computation of the 
partials of the acceleration due to the sun with respect to instantaneous posi- 
tion and velocity. 

— an executive routine calling various subprograms which evaluate the ac- 
celerations of a satellite due to the various forces acting on it. Modifica- 
tions were made to allow for the case of a two-body gravity model, and for 
the recomputation of the ephemeris quantities when necessary. 

— the subroutine DNVERT, which is a double precision matrix inversion 
routine using the Gauss -Jordan method of condensation with partial (column) 
pivoting, was modified to indicate an error condition due to a negative or 
zero pivot element. 
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4.2 New GEOSTAR-I ODP Modules and Significantly Modified NONAME ODP Modules 


The following GEOSTAR-I modules were written to allow for geopotential estimation, satisfy 
I/O interface requirements with the LUNGFISH MERGE and SOLVE programs, and numerically 
integrate the equations of motion and the variational equations. 

The orbit generator portion of the NONAME ODP program was replaced by a new set of subrou- 
tines which improves efficiency and provides additional capabilities to the system. Those NONAME 
subroutines which evaluate the force model associated with the satellite accelerations and accelera- 
tion partials have been left essentially unchanged, except for some additions which may improve 
the accuracy of the acceleration partials. 

Those subroutines effected by the new GEOSTAR-I orbit generator are summarized in Table 1. 


Table 1 


Replaced NONAME Orbit 
Generator Routines 

New GEOSTAR-I ODP Routines 
(approximate equivalents) 

COWELL 

INV2, CSTEP, MMATRX, 
SUMS, TEST,* CKDIFF* 

HERMIT 

HE MINT 

REARG, HHEMIT 

TABLEB* 

F, EGRAV, VEVAL 

FRCS, EGRAV, VEVAL 

INTGST 

TABLE, RK 

ORBIT 

ORBIT, EPHQAN, SWTEST 


*Variable stepsize version. 
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CSTEP 


Purpose : 

• To integrate the satellite equations of motion using a summed ordinate form of the 
Stormer/ Cowell predictor-corrector formulas for position, and Adams -Bashf or th/ 

Moulton formulas for velocity. 

• To integrate the variational equations using a summed "corrector-only" form of the 
Cowell formula for position partials and Moulton formula for velocity partials. 

Called By : 

ORBIT 

Calls: 

FRCS 

HEMINT 

VEVAL 

INV2 

MMATRX 

Method: 

• Integration of Equations of Motion 

Let t n+1 be the time point at which the satellite position and velocity is to be computed. Letting 
5L be the satellite acceleration vector at time t . , scaled by the factor h 2 , the predicted position 
and velocity vectors are computed from: 



i=0 


(Stormer) 



and successive corrections by: 


e 

3 + / CL* X 

n / i n+l~i 


i = 0 


(Bashforth) 


(Cowell) 
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X 


n + 1 


1 

h 




i = o 


a. m 


X n + 1 - i 


(Moulton) 


where x n+1 is computed using the subroutine FRCS whenever a successive correction is necessary. 
The maximum number of corrections allowed is three. After convergence, the sums I S n + 1 , 1 S n+ j 
are computed using: 


In the case when no drag is used, the velocity predictor is not used and the corrector is only ap- 
plied once. 

The nominal values for the stepsize h and order used for the integration of the equations of 
motion are : 

h = ICO sec 


Order - 11 (k e - 9) 


which can be changed on option. 

• Integration of Variational Equations 

Let t +1 be the time point at which the position and velocity partials are to be computed. Let 
y. , y . denote position and velocity partials with respect to a particular parameter at time, t i and 
y. the corresponding acceleration partial, scaled by h 2 . Note that this h 2 need not be the same as 
the one used in the satellite acceleration vector. 

The satellite position and velocity vectors at time t n+1 are first computed. If the stepsizes 
being used are not equal, these vectors are obtained by interpolation. Next, using the subroutine 
VEVAL, the acceleration partials with respect to position and velocity, 


<9x 


n + l 


dx 


and with respect to the parameter x., 


dx 


n + 1 


n + 1 ~ dXj 


are computed. 
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The vectors x n , v n are then formed by: 


• v 

h2 “O* f n + l + “Sn + 2Z + 


V 

n 


1 

h 


h2 4>* f„ 


+ IS n 



where I s n , 11 S n are the associated sums for the acceleration partials. The solution 



is then computed by 


y n + l 


[I ~ H] ’ 


X 

n 

v. 


( 6X6) 


where I is the 6x6 identity matrix and 

h2a 0* A n + l h2<X o‘ B n + 7| 

A n + 1 B n+1 

In the case when there is no drag, B n s o for all n, so that we simply have 

y n - [i -h 2 a* A n + 1 ] X n , 

and y n+1 computed directly from the corrector formula 


y n + l V n + h ^0 A n + 1. y n + l ' 

The above computations are performed for each position and velocity partial required. Note 
however, that the matrix (i -H)" 1 is independent of the particular partials being computed, so that 
it is only formed once. 

The nominal values for the stepsize and order used for the integration of the variational equation 

are: 

h = 100 sec 


Order - 7 (k y " 5) , 
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which can be changed on option, independently of those used for the integration of the equations of 
motion. 

• Integration Coefficients 

The integration coefficients are brought in through the COMMON/ AB CO E F/f r o m the BLOCK 
DATA subprogram. The Adams -Bashforth Predictor Coefficients (ALPHA), Adams-Moulton Cor- 
rector Coefficients (ALPHAS), Stormer Predictor Coefficients (BETA) and Cowell Corrector Coef- 
ficients (BETAS) for the summed ordinate forms of these, equations were computed to 18 digit ac- 
curacy (Reference 8). Two indices, IB(1) and IB(2), computed in the ORBIT subprogram indicate 
which set of coefficients to use for the equations of motion and which set of coefficients to use for 
the variational equations based on the requested or nominal order formula to be used in these 
integrations. 


Calling Sequence: 


CALL CSTEP (IEQ) 


COMMON Blocks Used: 


ABCOEF 

WORKER 

LIMITS 

XYZ 


Variables Not In COMMON: 


FORTRAN Name Format 

IEQ I 

SCO (3) D 

SPCO (3) D 

SN (6, 50) D 

C (6, 50) D 


ANPART 

GRBLOK 

DC 


Description 

Indicates which equations are to 
be integrated 

Nonsummed corrected value of 
position vector 

Nonsummed predicted value of 
position vector 

The storage array for the 
vectors X n and V n described in 
the method 

Instantaneous values of the 
position partials 
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EGRAV 


Purpose: 

• Calculates the satellite acceleration vector in rectangular coordinates due to the. geopotential. 

• Calculates acceleration partials in rectangular coordinates with respect to geopotential 
coefficients. 

Called By: 


FRCS 

VEVAL 


Method 


The accelerations due to the geopotential are given by: 


<9U _ <9U d g> 

dx (90 Sx 


where 0 = (r, x, 0) and u is the earth's potential. The components of the vector 


( 9u _ ( i9U ( 9U 

34, \ dl ' dK ’ 


3V 

(90 


are given by: 


(9U _ GM 
3 r r 2 


3 0 / a \ n n 

1 + 2^ (tJ 22 ( C n m COSmX + S nm sin mX ( n + 1 ). P n ( sin ^ 


rr-2 m = 0 


<9U 

dk 


30 , a v n n 

= ^ 22 w 22 ( Snm cos ,iA - Cnm sinm ^ tnp n ra < sin 


n = 2 m - 0 


(9U 

(90 


30 /a \ n n 

™ 22 (t) 22 cosmX + s ° m sinm ^ f p " +1 (sin ^ ■ mtan ^ p n ( sin 


n = 2 m=0 
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where 


GM = 
r = 
a = 

e 

K = 

0 = 

p m = 
r n 

c , s - 

n m ’ n m 


gravitational constant times mass of earth 
satellite radius vector 
earth's semi-major axis 
geocentric longitude of satellite 
geocentric latitude of satellite 
associated Legendre polynomials 
harmonic coefficients 


The matrix dcp/dx is given by: 


d± 

<?x 


by: 




~ x/r 

y/r 

z/r 


- y/s 2 

x/s 2 

0 

» S - 

xz/r 2 s 

- yz/r 2 s 

s/r 2 



s - (x 2 + y 2 ) 1/2 


To compute the Legendre polynomials P™ and the trigonometric terms, the following recursive 
relations are used: 


p n °(x) - [(2n-l) xP^.j (x)-(n*l)P n °_ 2 ]/n 
P: (x) = P™_ 2 + ( 2n - 1 ) (l-x 2 ) 1/2 P*Z 2 (x) , 


and for sectorials (m - n), since the first term is zero, 


P" (x) = ( 2n - 1) 00.. 


sinmA. 


cos mA. 


m tan i p 


2 cos A. sin (m - 1) A. - sin (m -2) A. 
2 cos A. cos (tn - 1) A. - cos (m - 2) A. 
(m - 1) tani/i + tan \p . 
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The gravitational accelerations in rectangular coordinates are then computed by: 


3U 

3x 


/au 

au 

au\ 

~<i£ 

Ur • 

a\ ■ 

d\jj) 

j9x_ 


The acceleration partials with respect to the harmonics C nm and S nm are obtained from 

9xJ 


and 


where 


3C 


a I 

'au 

au 

au’ 

= a c \ 

nm 

v a r - 

ax. ■ 

3i/» j 

(3 

au 

a 

au 

" Uc nra 

ar - 

ac 

nm 

3X. 


3 3 U 


3S„ 


)|] 


ra au a au a au\ [a^l 

= [as 3 r ' a s„ dK > <9 s U- 

\ nm nm nm 5 * 7 / ’X_] 


a 

au 

ac 

n m 

a r 


au 

ac 

n m 

a a. - 

a 

au 

ac 

nm 

dip 

a 

au 

as 

nm 

a r ~ 

d 

au 

as 

nm 

a>C - 

a 

au 

as 

nn 

3\p ~ 


GM / a e 

~ ( ” I (n + 1) cos m/VP™ ( sin \ji) 


- ~T~ ^ r ) nisinnAP” ( sin <p) 

GM / a e\ n 

“ ) cos mi. [P“ + 1 (sin.i//) - m tamp PJJ (sini/<)l 

5 3U GM W" 

as 3r - - 2 \ r / (n+ 1) smniVP” (smif) 

nm r * ' 


■ ? (y)” 


mcosmXP™ (sint//) 


GM 


~ sinmA. [P" + 1 (sin^) - mtan^PJ (sin</r)] 
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The Cand S coefficients are described by indices n andm. For the zonal harmonics, m = 0. Fo'r 
the sectorial harmonics, m = n. For the tesseral harmonics, m <n. Therefore, the C coefficients 
only fill one triangle of a matrix, as do the S. To conserve space, both sets of coefficients were 
combined into one matrix. A diagram of that matrix and the computations for the subscripts cor- 
responding to the n and m indices appear in Table 2 and a list of the SAO denormalized coefficients 
used by the program appears in the FMODEL COMMON block description in Section 5.1. 

Table 2 


Structure of SO x 33 
Harmonic Coefficients Array 






CS (30, 33) 




CO 


S!o° 

c29 
^3 0 


. S 1 Q 0 

& 3 0 b 3 0 


C« 

C 2 

C '. 

s 2 2 i 


S 29 


c 0 

29 



• . '■ t 


Sl 2 S 2 






n 2, 0 

S l SO 


C 3 0 

L 30 


C 2 9 

C 3 0 




Index 


Matrix Subscript Computation 

C coefficients 



n 



N 




m 



M + 1 

S coefficients 



n 



31 - N 




m 



33 - M 


Calling Sequence: 


CALL EGRAV (THETG, AE, GM, RASAT, DX, SWITCH) 


COMMON Blocks Used: 


FMODEL ESTGP 

XYZ SETSW 

GRBLOK IORBIT 

LIMITS 
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Variables Not in COMMON: 


FORTRAN Name Format 


THETG 

D 

AE 

D 

GM 

D 

RASAT 

D 

DX (3) 

D 

SWITCH 

L 


LAMBDA 

D 

P (32, 30) 

R 

C (30, 33) 

R 

S (30, 33) 

R 

PCNM ( 3)1 

D 

TO ! 

> D 

T1 

1 D 


j 


Description 

Right ascension of the Greenwich 
meridian 

Equatorial radius of earth 

Gravitational constant times mass 
of earth 

Right ascension of the satellite 

Satellite acceleration due to 
earth's gravity in rectangular 
coordinates 

A switch which if set to TRUE 
indicates that the acceleration 
partials with respect to the geo- 
potential coefficients will be 
computed 

X. -geocentric longitude of the 
satellite 

P n m -associated Legendre polynomials 
of degree m and order n 

c nm -coefficients of the cosine 
spherical harmonic 

S nra -coefficients of the sine 
spherical harmonic 

Intermediate quantities used in 
the computation of acceleration 
partials with respect to harmonic 
coefficients 
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EPHQAN 


Purpose: 

To control and determine when to calculate the three components of the moon's inertial unit 
vector and geocentric distance, the sun's inertial unit vector and geocentric distance, and the 
equation of the equinoxes. 


Called By : 

ORBIT 

Calls: 

COEFF 

DJUL 

EQN 

MOONAD 

SUN 

TDIF 

Method: 


In order to calculate the ephemeris quantities, the coefficients of nine 4th order degree poly- 
nomials are estimated using 11 calculated values of each quantity spaced over a 2.5 day arc as ob- 
servations. The degree of the polynomial, the number of observations, and the length of arc can be 
changed by altering these quantities in the BLOCK DATA subprogram. 


Calling Sequence : 


CALL EPHQAN 


COMMON Blocks Used: 


CONST 3 
COFIT 
COSAVE 
P COFIT 
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Variable s Not in COMMON : 


FORTRAN Name 

Format 

Description 

DJ 

D 

Julian date of ephemeris times 
used to compute the ephemeris 
polynomials 

AO (9, 11) 

D 

Array containing all ephemeris 
quantities to be used for the 
polynomials fit 

TIME (11)_ 

D 

Times of the ephemeris quantities 

VAR (11) 

D 

Ephemeris quantities to be used 
for the polynomial fit 

SI 

D 

S3 /2 

E 

D 

True obliquity of date 

CE 

D 

Cosine of true obliquity of date 

SE 

D 

Sine of true obliquity of date 
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GTIMIN/ GTIMOT 


Purpose : 

Measures the elapsed time of each iteration and computes the total elapsed time of a run in 
hundredths of seconds. 

Called By: 


MAIN 

Method: 

• The subroutine GTIMIN saves the time at which it is called but returns nothing to the pro- 
gram. However this subroutine must be called prior to the calling of GTIMOT. 

• The subroutine GTIMOT returns the elapsed time as measured by subtracting the time of 
day GTIMIN was called from the time of day GTIMOT was called. 

Calling Sequence: 


CALL GTIMOT (TX) 

TX— Name of the single precision variable containing the value of the elapsed time in seconds. 
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HEMINT 


Purpose : 

To interpolate for position, velocity and the associated partials. 


Called By: 


ORBIT 

CSTEP 


Method: 

The position and velocity components of the equations of motion and the associated partials 
are interpolated using the Hermite formula: 


where 


k 

f(x) = ^ h j (x ) f j + h j ^ i 

j=0 


h j ( x) - 

1 " 2 ( x " x i) ( 

x.) 4. 2 (x) 

h. (x) = 

(x- Xj ) ^(x) 



l. (x) - 

itj 

k J 

V (x) = l \ (x) 2] ' 

i=0 

i^j 


K 

n 


( x i - 


The order k is set to max (5, p/2). 


Calling Sequence : 


CALL HEMINT (TI, STEPZ, Kl, IEQ) 
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COMMON Blocks Used: 


WORKER 

LIMITS 


Variables Not in COMMON: 


FORTRAN Name 

Format 

Description 

TI 

D 

Interpolation time 

STEPZ 

D 

Distance between data points 

K1 

I 

Starting index of points to be 
interpolated 

IEQ 

I 

Indicates which equations are 
being interpolated 

M 

I 

Interpolation order (>5) 

SMI 

D 

Legendre interpolating polynomial 

SJMI 

D 

Derivative of Legendre interpolating 
polynomial 

HX 

D 

Hermite coefficient 

HXB 

D 

Derivative of Hermite coefficient 
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INV2/INV3 


Purpose/Method: 

To invert a matrix by using the Gauss-Jordan method of condensation with partial (column) 
pivoting. 

Note : 

The INV2 and INV3 subprograms are the same except that one passes variables through 
COMMON whereas the other passes variables using a calling sequence. The INV2 subprogram is 
used specifically to invert the matrix containing the Jacobian of accelerations with respect to posi- 
tion and velocity; the INV3 subprogram is used to invert matrices as required by other subprograms. 

INV2 Called By: 


CSTEP 


INV3 Called By: 


MAIN 

COfiFP 


Calling Sequence : 


CALL INV3 (N, A, M, DETERM) 


Variables: 


FORTRAN Name 

Format 

Description 

A (50, 50) 

D 

Matrix to be inverted and subsequent 
inverted matrix 

N 

I 

Dimension of portion of matrix to 
be inverted 

M 

I 

Maximum dimension of matrix to 
be inverted 

DETERM 

D 

Determinant of matrix inverted 
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MMATRX 


Purpose : 


To multiply an n x n matrix with an n x m matrix where n need not equal m . 


Called By: 


CSTEP 


Method: 


Standard matrix multiplication, 


C = A x B - 


where 




Cii = 71 3ik bk * 

i = 1 , 2, • • • n 
j = 1 , 2, • • ■ m . 

Calling Sequence: 

k = 0 

CALL MMATRX (A, 

B, N, M, C) 

Variables: 

FORTRAN Name 

Format 

Description 

A (6, 6) ] 
B (6, 50) j 

! 

D 

D 

Dimensions of input and output 
matrices. 

N, M 


I 

Input and output matrices 

C (6, 50) 


D 

Output matrix 
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ORBIT 


Purpose : 

To control the integration of the orbit and the variational equations. 


Called By : 

MAIN 

ORB1 

Calif 


SWTEST 

EPHQAN 

FRCS 

VEVAL 

1'ABLE 

SUMS 

CSTEP 

HEMINT 

REFCOR 

Method : 

Orbit has five functions: 1) to initialize the required variables and constants, 2) to compute 
the starting table for the integration, 3) to compute the initial first and second sums for the integra- 
tion, 4) to control the calculation of the ephemeris quantities, and 5) to control the orbit and var- 
iational equations' integration and interpolation processes. 


Calling Sequences: 


CALL ORBIT (PXPXO) 


COMMON Blocks Used: 


CONS 11 

CONST3 

ANPART 

ESTGP 

WORKER 

IORBIT 

NON2 

LIMITS 

CONVRG 

XYZOUT 

GRBLOK 

ABCOEF 

COFIT 
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Variables Not in COMMON: 


FORTRAN Name 

Format 

Description 

PXPXO (50, 6) 

D 

An array containing partials of po- 
sition and velocity with respect to 
parameters to be estimated 

BO 

R 

1/2 area of satellite/mass of 
satellite 

TI 

D 

Observation time in days 

IEND 

I 

Total number of array points for 
position, velocity, and position 
partials 

NO 

I 

Number of array points to reset 

IEQ 

I 

Equation indicator 

NEO 

I 


ETIME 

D 

Time to recompute the ephemeris 
quantities 

TOUT 

D 

Time from epoch to current in- 
tegration time in minutes 
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READGP 


Purpose : 

« To read in data to modify the stored values of the geopotential coefficients. 

• To read data determining which geopotential coefficients are to be estimated and to set up 
a table of the requested parameters. 

• To read other data options defining the geopotential estimation problem. 

Called By: 

OPT CARD 


Method : 

The call to READGP is initiated by the COEFGP option card. The remaining geopotential coef- 
ficient option cards are in two categories: EST and CHG cards. 

CHG cards contain the coefficient name and value to replace the prestored data base in COMMON/ 
FMODEL/CS. No other action takes place. There is no limit on the number of change cards. 

EST cards request the gravity coefficient indicated to be estimated. The initial estimate is 
taken from COMMON/FMODEL/CS unless the EST1 card contains a value for this estimate. The 
limit of EST cards is 50. The arrays of coefficient subscript indicators, estimated values and 
sigmas are sorted to place all C's first, then S's. Note that the a priori values and estimate values 
are initially the same. 

Calling Sequence: 

CALL READGP (INTP, LARRAY, DTI, T2, T3, T4, T5) 

COMMON Blocks Used: 


FMODEL 

ESTGP 

NON2 


Variables Not in COMMON : 

FORTRAN Name Format 

INTP I 

IARRAY (1) I 
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Description 

Data set to read from, usually FT05 

Resets logical indicator of ESTGP/ 
ISTATE 



IARRAY (2), (3) 

I 

Resets value of /ESTGP/ITERGP 

IARRAY (4) 

I 

Resets value of/NONl/JTEMP 

DTI 

D 

Not used 

T2 

D 

Sets logical indicator of/LIMITS/ 
ISWT (13) 

T3-T5 

D 

Not used 
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RK 


Purpose : 

To integrate the equations of motion and the variational equations to obtain the starting values 
of position, velocity, andpartials for the multistep Cowell integrator. 

Called By: 


TABLE 


Calls: 


FRCS 

VEVAL 


Method : 

The method used to obtain the starting values for the multistep integration of the equations of 
motion and variational equations is an eight order Runge-Kutta integrator. This method is defined 
as follows: Given a system of first order differential equations 

X = f(t, x) 


with initial vector x(t 0 ) - x 0 , the solution vector at time t x = t 0 +h, h a given stepsize, is com- 
puted by: 


x ( t i) = x ( fc o) + 840 [41 (K q +K 9 ) + 27(K 3 +K 5 ) + 272 K 4 + 216(K 6 +K3)J , 


where 


Kq ^ (*"0 ’ X o ) 

— I 4h _ 4 _ ^ 

K i - f ^0 + 27 - x o + W K o) 

K 2 = f(t 0 + TT.*o + llfo + 3Kj) 
K 3 = f(t 0 +X. *0 + 12 [K 0 + 3K 2 ]) 
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K 4 - f(t 0 +y. x 0 + ^ [K 0 + 3 K 3 ] ) 

K 5 = f(t 0 + T- *o + A [13K 0 -27K 2+ 42K 3 + 8K 4 ]) 

* 6 = f(to + lf’ *o + 4^0 [389 K 0 -54K 2 + 966 K 3 - 824K 4 + 243K S ]) 

K 7 = T^t 0 +h, x 0 + ^ [-231 K 0 + 81K 2 -1164 K 3 + 656K 4 - 122K S + 800K 6 ]j 
_ ( 5h 1 r 

K s = f + T ' *o + 288 L 127 K o + 18 K 2 “678 K 3 

+ 456 K 4 - 9 K s + 576 K fi + 4K ? fj 
K 9 = 7ft 0 +h, x 0 + gif) [1481 K 0 -81 K 2 +7104 K 3 

- 3376 K 4 + 72K s - 5040 K fi - 60 K ? + 720 K„ ] ) 

The stepsize h used in this integrator is nominally set at 24 seconds and can be modified on option. 
The subroutine accepts as input the values of the equations being integrated at time t . , the request 
time t. +1 = t, + h Cowell . It then performs N integration steps with the Runge-Kutta stepsize h, 
where N is the integer 


N 



followed by one final integration step with h = t. +1 -(t. +Nh). 


Calling Sequence: 


CALL RK 


COMMON Blocks Used: 


WORKER 

RKST 

RKT 

ANPART 

GRBLOK 
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Variables Not in COMMON: 


FORTRAN Name 

Format 

Description 

R* (10, 6, 50) 

D 

The Kj required by the method 

F (6, 50) 

D 

The f ; required by the method 


*NOTE: R is equivalenced to MYMD of COMMON/PRIORI/, thus writing over the flux tables after initial use. 
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SQLVGP 


Purpose : 

To solve the normal equations for the parameter corrections. This subroutine performs an 
accuracy check on the inverted normal matrix. If it is found that numerical difficulties were en- 
countered by the normal inversion process due to poor conditioning, a gradient method is used to 
obtain the parameter corrections. 

Called By: 


MAIN 


Calls: 


SYMMET 

DNVRT1 


Method : 

As input, SOLVGP is given the normal equations 

S Ax = b 


to solve for the parameter corrections Ax, where 

S = [ A T WA + P 0 “ 1 ] . 

the normal matrix summed over all measurements and 

b = A 1 ' WAO - P 0 -1 2 Ax, 


the right hand side of the normal equations. In these equations: 

A = the matrix of measurement partials with respect to the parameters to be estimated 
W = measurement weight matrix 
P 0 = a priori covariance matrix 
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AO = the vector of residuals 

SAx - the accumulated change in Ax over previous iterations. 

Reference the subroutine ESTIM (Reference 1, Vol. 2) for details on the formation of the normal 
equations. 

SOLVGP performs the following sequence of computations on the normal matrix S: 

(i) The matrix S is first normalized to form the matrix H where the elements of H are de- 
fined by: 



If s.j <0 for any i, transfer to step viii. 

(ii) Invert H using the subroutine DNVRT1; error exit from DNVERT resulting from zero di- 
visor defaults to the gradient option (viii). DNVRT1 uses the Gauss- Jordon method of con- 
densation with partial pivoting. The result of this inversion process is D where, hopefully, 
H' 1 = D. 

(iii) Test for the validity of the inverted matrix, H _1 ~ D by multiplying the matrix by its sup- 
posed inverse and comparing the result to an identity matrix i.e., compute the maximum 
element of the product matrix HD where 1 is subtracted from the diagonal elements. 

(iv) The resulting maximum, MAX, is tested against some tolerance, TOL, (initially 1 x 10 ~ 8 ). 
If MAX > TOL, then the resulting inverted matrix D is not accurately H" 1 and is con- 
sidered ill-conditioned. The gradient option (viii) is the default for this condition. 

(v) For MAX < TOL, the solution is valid and the elements of s _1 are given by 



(vi) Then the correction vector Ax is computed by Ax = s -1 b and S' 1 is the variance- covariance 
matrix. 

(vii) Program exit. 

(viii) Gradient option: This technique is an attempt to get a solution if errors resulted at (i), 

(ii) or (iv) above. 

For approximating the correction vector Ax, evaluate y = Sb, compute scale factor 

b T y 

cr = — ~ — » 
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and calculate x as x 


The error signals of 1, 2, and 3 are used for output to indicate to the user that the Ax correc 
tions obtained at that iteration are approximations resulting from the use of the gradient option. 
For error 3 the D matrix is available for display. 

Calling Sequence: 

CALL SOLVGP (DELTA, NPARAM) 

COMMON Blocks Used : 

PRIORI 

BEQ 

Variables Not in COMMO N: 

FORTRAN Name Format Description 

DELTA (50) D Parameter corrections 

NERR I Error indicator; zero = no error 
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STORGP 


Pur pose/ Method : 

The purpose of STORGP is to initialize the internal arrays of geopotential coefficients (CNM 
and SNM) used by the ORBIT subroutine VEVAL with the values from the data base reference. 
/FMODEL/CS. STORGP is called at the beginning of each DC iteration. Thus, in a DC geopotential 
estimation problem or when altering the geopotential model by the CHG car’d option, consistent 
values of the model are used. The arrays CNM and SNM contain values for the tesserals through 
order 3 and zonals through order 4. 

Called By: 

MAIN 


Calling Sequence : 


CALL STORGP 


COMMON Blocks Used: 


FMODEL 

CSVEVL 
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SUMS 


Purpose : 

To compute the first and second sums, r S n , n S n where n = order - 1, necessary for the summed 
form of the predictor -corrector formulas used in the multistep integration process. This is done 
only once for each iteration— during the initialization of the orbit generator. 

Called By 

ORBIT 


Method : 

Let denote the acceleration vector, or acceleration partial with respect to a parameter, 
where it has been scaled by the factor h 2 . The first and second sums are defined by: 

r S = V" 1 

n n 


V 1 *S = V' 2 * 


and are computed as follows. First, the sums n S n _ 1 and II S n _ 1 are computed by inverting the cor- 
rector formulas for x„ and x„ : 

n n 


*S _. - h x - y ' j8* x 

n~l n / , “ i r 

; = n 


k 



i = 0 


where a.* , /3.* are the Cowell and Moulton coefficients respectively. The required sums are then 
computed by 


r s 

n 

II 

I S . + x 

n - 1 

n s n 

II 

n s n _, + x s n 


Calling Sequence : 


CALL SUMS (IEQ, J3, J4) 
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COMMON Blocks Used: 


ABCOEF 

WORKER 

LIMITS 


Variables Not in COMMON : 


FORTRAN Name 

Format 

Description 

IEQ 

I 

Indicates which equation is being 



used 

J3 

I 

Indicates first equation to be used 

J4 

I 

Indicates last equation to be used 
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SUMTOB 


Purpose/Method : 

Moves the (K + I)-th parameter and its associated row of the normal equations: 

(SUMl) X (PARBUF) = (SUM2), with (LBLBUF) to the (L + J)-th position of the normal equa- 
tions in the B matrix form 

(BMATRX) X (PARAM) = (BRHS), with (LABEL). 


Called By: 


WTBMAT 


Calling Sequence: 


CALL SUMTOB (N, K, I, L, J) 


COMMON Blocks Used: 


PRIORI 

BEQ 


Variables Not in COMMON: 

FORTRAN Name Dim. Format 

N II 

K 1 I 

I 1 I 

L II 

J 1 I 


Description 

Total number of parameters (=matrix 
size) 

Number of parameter of all preceding 
parameter groups in the input block SUM2 
(0 <K< N) 

Position of the parameter s of the 
current parameter group to be moved 
(1 < K + I < N) 

Number of parameters in output param- 
eter block PARAM (0 < L < N ) 

Position of the parameter in the cur- 
rent parameter group in B-matrix 
order (1 < L + j < n) 
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TABLE 


Purpose: 

To form the initial table of starting values for the equations of motion and the variational 
equations to be used in the multistep Cowell integrator. 

Called By: 


ORBIT 


Calls: 


RK 


Method: 

The computation of the starting values required by the multistep integrator is done by using a 
single step 8 th order Runge-Kutta integration method. The algorithm used to form the starting 
tables is the following: 

Case I: Equations of motion and variational equations being integrated with the same stepsize. 
In this case, the total number of time points produced for the starting tables is 

K = max (Nj, N 2 ) , 


where 


= ?! - 2 ; 

Pj and P 2 the orders of the integrators to be used for the equations of motion and variational equa- 
tions respectively. The required times t. = t 0 + ih, i = 1, 2, • • • k are computed and used as in- 
put to the subroutine RK together with the initial values to compute the necessary starting table. 

Case II: Equations of motion and variational equations being integrated with different stepsizes. 

In this case, the total number of time points produced for the starting tables is 

K = Nj + N 2 . 

where n. is defined as above. If h x and h 2 are the stepsizes to be used to integrate the equations of 
motion and the variational equations respectively, then the required times are given by 
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+ ih 1 

i = 1, • 


+ jh 2 

II 

I-* 

’ ' n 2 


and are used as input to the subroutine RK together with the initial values to obtain arrays con- 
taining the required starting tables. The starting tables required at the two stepsizes are then ob- 
tained by sorting the values contained in the RK produced arrays. 

Because of core storage limitations, the restriction 

p i + p 2 < 22 

is required whenever two different stepsizes are used. 

Calling Sequence: 


CALL TABLE 


COMMON Blocks Used: 


WORKER 

RKT 


Variables Not in COMMON: 


FORTRAN Name 
IND (30) 


K 

RT1 


RT2 


Format 

I 


I 

I 


I 


Description 

Array of indicators to designate 
which starting points belong to the 
equations of motion and which 
starting points belong to the 
variational equations 

Total number of starting points 
needed 

Time at which starting points are 
to be computed for the equations 
of motion 

Time at which starting points are 
to be computed for the variational 
equations 
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IRTl 


IRT2 


ITR 


I Number of starting points to be 

computed for the equations of 
motion 

I Number of starting points to be 

computed for the variational 
equations 

I Indicates which equations are being 

re-organized 
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SWTEST 


Purpose : 

To set internal switches which: 

• identify which perturbative effects are being applied to the satellite accelerations and 
which perturbation effects are being applied to the computation of the position partials; 

• indicate what perturbations have to be recomputed if the variational equations and the 
equations of motion are being integrated with two different stepsizes; 

• indicate which variational equations are being integrated. 

Called By: 


ORBIT 


Method : 

The routine defines an array of indicators which are used to define computed GO TO statements 
in the force computations of the acceleration partials. The values for the indicators are determined 
by: 


(i) whether the equations of motion and variational equations are being integrated together 
or separately; 

(ii) what forces are to be applied to the variational equations; 

(iii) what forces are to be applied to the equations of motion. 

Calling Sequence: 


CALL SWTEST (SWITCH) 


COMMON Blocks Used: 


WORKER 

ESTGP 

LIMITS 

SETSW 
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Variables Not in COMMON: 


FORTRAN Name 
SWITCH 


Forma t Description 

L A switch which if set to TRUE in- 

dicates that the equations of 
motion and the variational equations 
are being integrated with the same 
stepsize; if the switch is set to 
FALSE, it indicates that the equa- 
tions of motion and the variational 
equations are being integrated with 
two different stepsizes. 
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VEVAL 


Purpose : 

To compute the partials of the satellite acceleration vectors with respect to the position and 
velocity and, if present, emissivity, drag, and harmonic coefficients. 


Called By : 

ORBIT 

RK 

CSTEP 

Calls : 

REFCOR 

EGRAV 

DENSTY 

DOTPRD 

Method : 

The acceleration partials are computed from accelerations of the form 


— 4- F . 4- F 4- F + F 

g— drag sr sun moon 


where 

u = the earth's geopotential 
F drog = acceleration vector due to atmospheric drag 
F sr = acceleration vector due to solar radiation pressure 
F sun = acceleration vector due to solar gravity 
F mo on = acceleration vector due to lunar gravity 

In the current NONAME ODP version of this subroutine, the acceleration partials include ef- 
fects of the earths potential to fourth order terms, drag and lunar gravity. The details of these 
computations can be found in Reference 1, Vol. 2. The GEOSTAR- 1 modification of this subroutine 
can be summarized as follows: 
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(i) effects due to solar gravity and radiation pressure have been included in the computation 
of acceleration partials. These solar gravity effects are given by: 




dx 

^y 



dF 

sun _ 
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where 

GM s = gravitational constant times mass of sun in earth masses, 
x s = (x s , y s , 2 S ) = solar position vector, 

and the solar radiation effects are given by: 


"dx 

s r 

dx 

s r 

dx 

s r 

dx 

3y 

d Z 

dy'sr 

dy 

J s r 

dy 

J s r 

dx 

d y 

dz 

dz 
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dz 
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dx 

dy 

dz 
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where 

C R = satellite reflectivity constant 

A s satellite cross-sectional area 

M = satellite mass 

(ii) Computations have been rearranged so that any of the particular effects acting on the ac- 
celeration partials can be excluded on option. Only an option to compute two-body partials 
has been included. 

(iii) The subroutine computes partials of acceleration with respect to harmonic coefficients 
by calling the EGRAV subroutine. The formulation is documented in that subroutine. 

(iv) Computations have been arranged so that if the variational equations are being integrated 
with a stepsize different from that used for the equations of motion, the intermediate 
quantities which are used in the computation of the acceleration partials (normally brought 
in through COMMON) are recomputed in VEVAL. 

(v) To allow for possible variations in the low order harmonics currently used in the ac- 
celeration partials computation, these coefficients are now initialized by the STORGP 
subroutine and brought into VEVAL through COMMON (CSVEVL). 
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Calling Sequence : 


COMMON Blocks Used : 

CONST1 

CONST3 

DRGBLK 

CSVEVL 

SETSW 

Variables Not in COMMON : 

FORTRAN Name 
RASAT 
THETG 

LAMBDA 

ELEM(6) 

VARD(9) 

EQ 

CONS 

DPXUVM 

RRSUN 

DPXUV 


CALL VEVAL 


VRBLOK 

FLX'BLK 

ESTGP 

GRBLOK 

MOGNGR 


LIMITS 

COFIT 

XYZ 

anpart 


Format 

D 

D 

D 

D 


D 


D 


Description 

Right ascension of the satellite 

Right ascension of Greenwich at 
integration time 

Satellite longitude 

Satellite position and velocity- 
vectors referenced to the true equator 
and equinox of date 

An array containing lunar and solar 
position vectors and vector 
magnitudes 

Equation of equinox 


R Constant used in the computation 

of the acceleration partials with 
respect to the drag parameter 

D Dot product of satellite and lunar 

position vector 

D Square of lunar or solar position 

vector 

D Dot product of satellite and solar 

position vector 
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P(3), PMAG(3) D 

SUNMAT (3,3) D 

SOLMAT (3,3) D 

SLUMAT (3, 3) D 


Vector and vector magnitude used 
in the sunlight determination 
computation 

Matrix of the satellite accelera- 
tion partials due to solar gravity with 
with respect to instantaneous 
position 

Matrix of the satellite accelera- 
partials due to solar radiation 
with respect to instantaneous 
position 

Matrix of the satellite accelera- 
tion partials due to lunar gravity 
with respect to instantaneous 
position 
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WTBMAT 


Purpose 

To output on the specified device the normal equations, together with the latest values and the 
labels of the parameters, in the B matrix format. 

Called By : 

MAIN 


Calls : 

SUMTOB 

SYMMET 


Method : 

The parameter sequence in the normal equations in GEOSTAR-I ODP is different from that 
required by the B matrix format. Figure 2, Section 3.3.1 shows this parameter order. WTBMAT 
takes each parameter type group and assigns labels and sorts the parameters into the B matrix 
format within that particular group. 

The subprogram SUMTOB is used to shift rows of a particular parameter type group into the 
designated B matrix areas. After all parameter groups are sorted and labeled, the normal equa- 
tions matrix is written by rows in the B matrix format. 

The total variance, as required by the SOLVE program, is computed as: 

VI = (RMS) 2 * [NOB - (NPARAM - 1)] , 


where 


RMS 

NOB 

NPARAM 


total weighted standard deviation 

total number of observations used 

total number of parameters being estimated. 


Calling Sequence: 
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CALL WTBMAT (RMSTOT, NOWTOB, NSTEST, NBIAS, 
NPARAM, BNAME, IDMAT, BMT) 



COMMON Blocks Used: 


CONST 1 STANUM 

ESTGP PRIORI 

CELEM BEQ 


Variables: 


FORTRAN Name 

Format 

Description 

NOWTOB 

I 

NOB, number of weighted observa- 
tions used in computing RMSTOT 

NSTEST 

I 

Number of stations whose locations 
are to be estimated 

NBIAS 

I 

Number of instrumental bias 
parameters 

NPARAM 

I 

Total number of parameters to be 
estimated 

BNAME (3) 

R 

B matrix name (12 EBCDIC 
characters) 

IDMAT 

I 

B matrix identification number 
' (<99999) 

BMT 

I 

Output data set reference number 

VI 

D 

The total variance 

RMSTOT 

R 

Total weighted standard deviation 
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4.3 Modifications to Existing SOLVE, EIGENVALUE Modules 

The modifications made to the SOLVE subroutines are designed to: 

• Call the new subroutines. 

• Allow for 4-digit station code labeling. 

• Change the punched card formats to conform to the GEOSTAR-I ODP input requirements. 

• Accept station location parameters in rectangular coordinates and convert them to spheri- 
cal coordinates for printing and punching. 

A summary of these modifications for the GEOSTAR-I SOLVE follows: 

• CALTYP— Modified to recognize 4-digit station numbers, and to adjust its label count 
accordingly. 

• INVERT— Modified to call for the pseudoinverse solution (ANDREE) on option, with cor- 
responding changes in matrix handling, etc. 

• LBLSUP— Modified to recognize 4-digit station numbers and their special suppression 
features, and to zero out the proper labels accordingly. 

• O PARC— Modified to punch out updated state parameters in a format suitable for input back 
into ODP for a GEOSTAR-I iteration. 

• OPGRAV— Modified to punch out updated harmonics in a format suitable for input back into 
NONAME for a GEOSTAR-I iteration. 

• OPSTAT— Modified to recognize 4-digit station parameters; accept B matrix input given in 
rectangular coordinates; output rectangular or spherical coordinates as required; punch 

out STAPOS cards in a format suitable for input back into the ODP for a GEOSTAR-I iteration. 

• SUPRSS— Modified to recognize an input CB (combined) matrix. 

• UPCOMB— Modified to recognize 4-digit station numbers. 

In addition the SETEIG subroutine in the EIGENVALUE program has been modified to call the 
subroutine ELIM which eliminates the station position parameters from the input matrix. 
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4.4 New GEOSTAR- 1 SOLVE Modules and Significantly Modified LUNGFISH SOLVE Modules 

The following GEOSTAR- 1 modules were written to allow for pseudoinversion, to satisfy the I/O 
interface requirements with the GEOSTAR-I ODP program and to increase the "combined” matrix 
handling capabilities of the original SOLVE program. The subroutine MAIN, the SOLVE control 
program, was modified considerably to achieve the above requirements and hence these modifica- 
tions are detailed in the following pages, as well as the new modules. 

The modifications to MAIN are now detailed. 

The first section of MAIN reads data cards one through five. In this section, modifica- 
tions have been made as follows: 

(i) Three additional options, IOPT5, IOPT6, and IOPT7 have been included in the read-in of 
the option data card, data card #2. These options are: to allow the input of a combined 
matrix with associated backsubstitution and parameter set matrices (IOPT5); to compute 
the Penrose pseudoinverse instead of a normal inverse (IQPT6; the ANDREE algorithm is 
used); and to accept incoming B matrices as representing the same arc (although different 
data types), allowing for the combining of the entire matrices, including state parameters 
(IOPT7). 

(ii) In accordance with IOPT5, input tapes 19 or 29 (containing input parameter sets and 
input backsubstitution matrices to go with an input combined matrix) are repositioned 
so that they may be used in subsequent parts of the program. 

The second section of MAIN reads data cards six and seven, B matrix information, and then 
finds and transfers the B matrix from input tape 18 to work tape 28 (suppressing the B matrix as 
required). 

The following changes have been made in this section: 

(i) Provision has been made for reading in optional values for AED and FLAT (parameters 
used in the conversion of rectangular to geodetic coordinates). These values are read in 
on B matrix data card 5; if the appropriate fields are left blank, nominal values for AED 
and FLAT are used. 

(ii) For the combine arcs option (IOPT7), provision was made for identifying such matrices 
by setting their ITYPE value in the matrix header to 1, thus identifying them as B 
matrices. The combined arcs matrix, when finally output, will be for all functional 
purposes equivalent to a B matrix. 

(iii ) Provision was made to recognize an input combined matrix (IOPT5). 

(iv) Coding was added to write a negative dummy record at the end of unit 19 once all the 
parameter sets were read in. This is for use in OPSTAT, OPGRAV and OPARC. 
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The third section of MAIN deals with the elimination of arc parameters from the input matrices. 
The following change has been made: 

If this is a combine arcs run (IOPT7), then switches are set so as not to eliminate any arc 
parameters. Thus the arc parameters will be combined along with the rest of the 
matrix. 

The fourth section of MAIN deals with the combining of the various arc-eliminated B matrices. 
The following changes have been made: 

(i) Coding from ELIM has been placed in MAIN to write out a negative dummy record only 
once on unit 29, at the end of all the backsubstitution. 

(ii) Coding has been added to complete implementation of the combined arcs option (IOPT7). 
Once the combining of matrices has been completed, then, a parameter set is read in from 
tape 19 and added immediately after the combined matrix on unit 28. The resulting matrix 
on tape 28 is then a B matrix, suitable for reintroduction into SOLVE, input to MERGE, etc. 

(iii) Depending on the exact value of IOPT7, the option has been introduced either to stop at 
this point and output the combined arcs matrix on tape 28, or to continue on to invert and 
print out final solutions. 

The fifth section of MAIN inverts the matrix by calling INVERT and its associated subroutines. 

The following change has been made in this part of MAIN: 

The call statement for INVERT now includes IOPT6 for use in INVERT (INVERT calls 
ANDREE for the pseudoinverse solution if IOPT6 is on). 

The sixth section of MAIN deals with the output. 

The following change has been made in this section: 

A mechanism was introduced for handling standard and optional values of AE and FLAT. 

The seventh and last major section of MAIN deals with backsubstitution and the printout of 
related output. 

The following change has been made in this section: 

If IOPT7 is non-zero, i.e. if this is a combined arcs run where no arcs were eliminated, all 
subsequent backsubstitution procedures are skipped. 

The very last sequence edits whatever matrices have been designated. One change has been 
added: 

A data card is read in; if it is not 77777 in value, then the program loops back to start, be- 
ginning with the first data card. 
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ANDREE 


Purpose : 

To compute the pseudoinverse of a given matrix. 
Called By: 


INVERT 


Method: 

The Andree algorithm is used to obtain the pseudoinverse from an input matrix. By 
pivoting about each largest diagonal element and discarding those diagonal elements which 
are unacceptably small (in the ill-conditioned case), an S-matrix is obtained so that SBS t = E, where 
B is the input matrix and E is a diagonal matrix of unitary or zero elements. If no elements were 
discarded in the pivot search, then E is the identity matrix and B" 1 is obtained by a normal inversion, 
S T S = B' 1 . 

If some elements were discarded however, then the algorithm computes a U matrix from the 
S matrix. Non-diagonal elements are set to the negative of their corresponding values in the 
S matrix if the diagonal for that row was eliminated, i.e., is zero in the E matrix. Otherwise, non- 
diagonal elements are zeroed, and the diagonal elements of the U matrix are equal to the diagonal 
elements of the E matrix. 

u T bu is formed and all rows and columns that are zero in the result are deleted. The re- 
mainder becomes the upper left portion of a matrix whose other portions are zero matrices: this 
is our new matrix C # , and B # , the pseudoinverse equals uc # u T . 

Calling Sequence : 

CALL ANDREE (V, N, NR, EPS, B, U, R, ISV) 


Variables : 



FORTRAN Name 

Format 

Description 

V(140, 140) 

D 

The matrix to be pseudoinverted 

N 

I 

The dimension of V 

EPS 

R 

Number of inaccurate digits, used in 
determining whether to zero a given 
pivot element 
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ISV 

I 

Work tape unit number (set at 25) 

NR 

I 

The computed rank of the matrix 

B(140) 

D 

Work array 

U(140, 140) 

D 

Work array 

R(140) 

D 

Work array 


The pseudoinverse of the input matrix is placed in the input array V for output. 
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OUTRAD 


Purpose : 

Converts radius to angular or time equivalents. 
Called By : 

OBSTAT 


Method : 

Proper combination of integer and real number arithmetic for the separation of higher integer 
units from remainders. 

Calling Sequence: 

CALL OUTRAD (RAD, IH, IM, S, K) 


Variables: 



FORTRAN Name 

Format 

Description 

RAD 

D 

Radian input 

IH 

I 

Equivalent number in hours or 
degrees 

IM 

I 

Minutes 

S 

D 

Seconds 

K 

I 

Input indicator specifying time or 
angular conversion 
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PHLINN 


Purpose : 

To convert tracking station geodetic rectangular coordinates to geodetic spherical coordinates. 
Called By : 

OPSTAT 

Calls: 

DARCTN 

Method : 

An iterative technique is used to compute geodetic latitude and elevation above spheroid from 
geocentric x, Y, Z of the input station. See Figure 3. 

Initially set 

t = e 2 Z . 

Then for each iteration solve for t as follows: 
z t = z + t 

N = (x 2 + Y 2 + Z t 2 ) 1/2 

sin cf> = Z t /N 
v = a/^1 - e 2 sin 2 

t = ve 2 sin<£ . 

When t converges, compute 



<f> - sin -1 4> (latitude) 


H - N _ v (height) 
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and 


X. = tan -1 (Y/X) 


Calling Sequence: 

CALL PHLINN (M, X, Y, Z, 


Variables 

FORTRAN Name Format 


M 

I 

X 


Y 

D 

Z 


AE 

D 

FLAT 

D 

PHI 

D 

XAMBDA 

D 

H 

D 


(longitude) . 

AE, FLAT, PHI, XAMBDA, H) 

Description 

Station identification indicator 

I Rectangular coordinates of 
1 station position 

Semi- major axis of reference 
ellipsoid 

Flattening coefficient of reference 
ellipsoid 

Geodetic latitude 

Geodetic east longitude 

Height above spheriod, in meters 
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V COMMON BLOCK VARIABLE DESCRIPTION 


The following sections contain the COMMON block descriptions of the COMMON areas used in 
the GEOSTAR-I ODP and SOLVE programs. These descriptions include the variables contained 
in these areas, their type and meaning, and the subroutines which either define or use these vari- 
ables, for those which are not used in a major number of the routines. Also, following each section 
of COMMON block variable descriptions, a section containing a cross reference table of all the 
COMMON areas as used in each subroutine is presented. 
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5.1 NONAMfi — GEOST AE.-I ODP COMMON Blocks 

This section contains a description of all the COMMON areas used in the GEOSTAR-I ODP 
program. 
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/ABCOEF/ 


COMMON /ABCOEF/ALPHA(102), ALPHAS (102), BETA(102), BETAS(102), IB(2) 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

ALPHA(102) 

D 

Adams -Bashforth predictor 
coefficients for orders 4-15 

BLOCK 

DATA 

CSTEP 

SUMS 

ALPHAS (102) 

D 

Adams-Moulton corrector 
coefficients for orders 4-15 

BLOCK 

DATA 

CSTEP 

SUMS 

BETA(102) 

D 

Stormer predictor coefficients 
for orders 4-15 

BLOCK 

DATA 

CSTEP 

SUMS 

BETAS (102) 

D 

Cowell corrector coefficients 
for orders 4-15 

BLOCK 

DATA 

CSTEP 

SUMS 

IB(2) 

I 

IB(1)— starting point for the set 
of integration coefficients to be 

ORBIT 

CSTEP 

SUMS 


used for the equations of motion 
integration 

IB(2)— starting point for the set 
of integration coefficients to be 
used for the variational equa- 
tions integration 
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/ANPARty 


COMMON / ANPART/ PRDD (3 , 6) 





Program 

Program 




Where 

Where 

Variable 

Type 

Description 

Defined 

Used 

PRDD(3, 6) 

D 

Matrix of the satellite accelera- 

VEVAL 

CSTEP 



tion partials with respect to 


ORBIT 



instantaneous position and 


RK 



velocity. The matrix is of the 


TABLEB 


form: 


dx 

dx 

dx. 

f 1L 

dx 

iiT 

dx 

3y 

dz 

| dx 

dy 

<9z| 

d 'y 

dy 

dy 

! a 

dy 

dy 1 

nr 

dy 

TT 

1 dk 

i 

dy 

dz J 

d'z 

d'i 

dz 

i 

i dz 

d'z 

d'z | 


W 

d z 

[_ dk 

dy 

d zj 


where those quantities within 
the dotted brackets are computed 
only if drag is applied. 
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/BEQ/ 


COMMON /BEQ/SUM2 (50), PARBUF{50), LBLBUF(50), BMATRX(50, 50), 
BRHS(50), PARAM<50), LABEL(50), IPARAM(50) 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

SUM2(50) 

D 

Right hand side of normal 
equations in GEOSTAR- 1 
parameter order 

ESTIM 

ESTIM 

WTBMAT 

SUMTOB 

SOLVGP 

PARBUF(50) 

D 

Buffer for parameter values 
in GEOSTAR- 1 order, by row 

WTBMAT 

WTBMAT 

SUMTOB 

LBLBUF(50) 

I 

Numerical identifiers for the 
type of parameters in PARBUF, 
in the same order 

WTBMAT 

WTBMAT 

SUMTOB 

BMATRX{50, 50) 

D 

Symmetrical matrix of the 
normal equations in B matrix 
form, for multi-arc processing 

WTBMAT 

WTBMAT 

SUMTOB 

BRHS(50) 

D 

The associated right hand 
side of BMATRX 

WTBMAT 

WTBMAT 

SUMTOB 

PARAM(50) 

D 

Parameter values in 
B matrix order by row 

WTBMAT 

WTBMAT 

LABEL(50) 

I 

Numerical identifiers for the 
type of parameters in PARAM, 
in the same order 

WTBMAT 

WTBMAT 

SUMTOB 

IPARAM(50) 

I 

Array used by subroutine 
SUMTOB to order the param- 
eter values 

SUMTOB 

SUMTOB 
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/CELEM/ 

COMMON /CELEM/ELEMST(6), ORBELA(6), XNV, EC 


Variable 


Type 


ELEMST(l) 

ELEMST(2) 

ELEMST(3) 

ELEMST(4) 

ELEMST(5) 

ELEMST(6) 




> 


j 


D 


ORBELA(l) 

ORBELA(2) 

ORBELA(3) 

ORBELA(4) 

ORBELA(5) 

ORBELA(6) 

XNU 

EC 


>1 


> 


D 


D 

D 


Description 

Orbital Elements - Inertial 
Position and Velocity Vectors 

X - meters 

Y - meters 
Z— meters 

X - meters/second 

Y - meter s/second 
Z - meters/second 

Orbital Elements - Osculating 
Keplerian 

Semi- major axis - meters 

Eccentricity 

Inclination - adians 

Longitude of ascending node - radians 

Argument of pericenter - radians 

Mean anomaly - radians 

True anomaly - radians 

Eccentric anomaly - radians 


Program 

Where 

Used 

ELEM 


POSVEL 


ESTIM 


Purpose 

Input - values of ELEMST(6) (inertial position 
and velocity vectors) from MAIN 
program 

Output - return values of ORBELA(6) (osculating 

orbital elements), XNV and EC to MAIN 
program 

Input - values of ORBELA(6) (osculating orbital 

elements) from MAIN program 

Output - values of ELEMST(6) (inertial position 

and velocity vectors) to MAIN program 

Input - values of ELEMST(6) (inertial position 

and velocity vectors). 
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/ CGEOS/ 


COMMON / CGEOS/ DAYREF, DAYSTP, DAYSTA, ISATID, SIG1, SIG2, 

MTYPE, NMEAS, ISTA 

NOTE 

Subprograms GEOSRD or DODSRD are utilized by the system to read records of ob- 
servational data, select and delete records on request, and preprocess optical and 
time observational data as requested by various indicator variables. Part of the final 
observational data is returned to MAIN through COMMON/CGEOS/. The remainder 
of the observational data is returned to MAIN through COMMON/ PREBLK/. Each 
record of observational data may contain one or two observations. 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

DAYREF 

D 

Reference data in days since Jan. 
0.0 in A1 time 

MAIN 


DAYSTP 

D 

Stop date for the DC arc in days 
since Jan. 0.0 in A1 time 

MAIN 


DAYSTA 

D 

Epoch of observations in 
days elapsed from Jan. 0.0 
UTC of the year of the 
epoch of the current 
run. 

GEOSRD 
DODSRD 
Converted to 
A1 time in 
MAIN 

main 

ISATID 

D 

Satellite identification code 

GEOSRD 

DODSRD 

MAIN 

SIG1 

R 

Standard deviation of 1st 
observation in this record 

GEOSRD 

DODSRD 

MAIN 

SIG2 

R 

Standard deviation of 2nd 
observation if the record 
contains two observations 

GESORD 

DODSRD 

MAIN 
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/CGEOS/ (continued) 


Var: 

MTYPE 


NMEAS 

ISN 


NOTE 

Units of SIG1 and SIG2 are : 


right ascension: 
declination: 
azimuth: 
elevation: 

X angle: 

Y angle: 
range: 
range rate: 
direction cosines: 


seconds of arc 
seconds of arc 
seconds of arc 
seconds of arc 
degrees 
degrees 
meters 
meters/sec 
dimensionless 


able 


Type 


Description 


Program 

Where 

Defined 


I 


Observation Type in this record 

= l-a. and s (right ascension 
and declination) 

= 2 - r (range) 

= 3 - r (range rate) 

= 4 - Doppler (converted to 
range rate) 

= 5 - i , m (Minitrack direction 
cosines) 

= 6 - X, Y angles 

= 7 - azimuth and elevation 


GEOSRD 

DODSRD 


I Number of observations in this 

record 

NMEAS = 1 for MTYPE = 2, 3, 4 

= 2 for MTYPE = 1, 5, 6, 7 

I GSFC code number for station 

recording observations in this 
record 


Program 

Where 

Used 


MAIN 

PREDCT 


MAIN 


MAIN 



/COFIT/ 


Variable 
C(5, 9) 


C(l, 1)-C(5, lfl 
C(l, 2)-C(5, 2) V 
C(l, 3)-C(5, 3) J 

C(l, 4)-C(5, 4) 

C(l, 5)-C(5, 5)' 
C(l, 6)-C(5, 6) V 
C(l, 7)-C(5, 7) J 

C(l, 8)-C(5, 8) 


C(l, 9)-C(5, 9) 


COMMON /COFIT/C (5, 9), DAY1, CENTER, H4 


Type 


Definition 


Program 

Where 

Defined 


Program 

Where 

Used 


Coefficients of a 4th degree 
polynomial fit to nine ephemeris 
quantities. These coefficients 
are estimated by using 11 values 
of each quantity calculated by 
subprograms MOONAD and 
SUN equally spaced over a 2.5 
day interval starting at DAYl. 


EPHQAN 


FRCS 

VEVAL 


Specifically C represents coefficients for: 


D 

X m 

EPHQAN 

FRCS 

VEVAL 

D 


EPHQAN 

FRCS 

VEVAL 

D 

x s 

EPHQAN 

FRCS 

VEVAL 

D 


EPHQAN 

FRCS 

VEVAL 

D 

E 

q 

EPHQAN 

FRCS 

VEVAL 

PREDCT 
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/COFIT/ (continued) 


Variable 


DAY1 


CENTER 


114 


Type 


Description 


Program 

Where 

Defined 


D Days elapsed since January ORBIT 

0.0 of the epoch year of the 
initial elements 

D = DAY1 + 1.25 days EPHQAN 


I Order of ephemeris polynomials EPHQAN 

where 


Program 

Where 

Used 


FRCS 

VEVAL 


FRCS 

VEVAL 

PREDCT 

FRCS 

VEVAL 


X = Unit vector to moon re- 

m 

ferred to true equator and 
equinox of date 

R = Distance to moon - meters 

m 

X = Unit vector to sun referred 

s 

to true equator and equinox 
of date 

R s = Distance to sun - meters 

E q = Equation of the equinoxes - 
radians 
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/ CONST 1/ 


COMMON /CONSTl/THETG0(10), THDOT1, GM, AE, FINV, CD, ASAT, 
MSAT, ADDR, SRAD, EMISS, MSUN, MMOON, RPRESS 


Variable 

Type 

Description 

THETGO(IO) 

D 

Greenwich Mean Sidereal Time in radians at 0 hours UTl on 
January 0 for the years 1960-1969 

THDOTl 

D 

Relative increase of angle between mean Greenwich meridian 
and the mean equinox of date - radians per 24 hours UTl time 

GM 

D 

GM e - meters 3 /seconds 2 

AE 

D 

R„ - meters 

FINV 

R 

1/f 

CD 

R 

c D 

ASAT 

R 

A s - meters 2 

MSAT 

R 

M s - kilograms 

ADDR 

T 

C D estimation indicator. If ADDR >0, C D will be estimated. 

SRAD 

I 

C R estimation indicator. If SRAD >0, C R will be estimated. 

EMISS 

R 

c R 

MSUN 

R 

M 0 / m e 

MMOON 

R 

M„/M e 

RPRESS 

R 

p 0 


where 

G = universal gravitational constant 
M E = mass of sun 
M H = mass of moon 
M s = mass of satellite 
A s = area of satellite 


f = flattening of earth's reference ellipsoid 

R E = semi-major axis of earth's reference ellipsoid 

C D = satellite drag coefficient 

C R = satellite reflectivity coefficient 

P Q = solar radiation pressure in the vicinity of the earth. 



/C0NST1/ (continued) 


THETGO(l) = 
( 2 ) = 

(3) = 

(4) = 

(5) = 

( 6 ) = 

(7) = 

( 8 ) = 
(9) = 

(10) = 

THDOT1 

RPRESS 


Value 

1.722 186 300 radians 

1.73 5 2 2 2 6 56 radians 

1.73 1 05 6 239 radians 

1.726 889 824 radians 

1.722 723 442 radians 

1.73 5 7 5 9 816 radians 

1.73 1 593 3 99 radians 

1.727 427 000 radians 

1.723 260 602 radians 

1.73 6 2 9 6 992 radians 

.017 202 791 266 radians per 24 

.45 x 10" 5 newtons/meter 2 


For the year 

1960 

1961 

1962 

1963 

1964 

1965 

1966 

1967 

1968 

1969 

UTl time 


The following nominal values may be changed by optional input cards : 


GM„ = 3.986 032 x 10 14 meters 3 / seconds 2 
R e = 6378165. meters 
1/f = 298.252 

C c = o 

A s = 0 
M = 0 

S 

C R = 0 

M q /M e = 332951.3 
M/M e = .0122999 


ADBR and SRAD axe used to request adjustments to C D and C R respectively, and are initially de- 
fined to be 0. If it is desired to adjust C D , option card DRAG is read to redefine ADDR to be +1. 
If it is desired to adjust C R , option card SOLARD is read to redefine SRAD to be +1. 


The CONSTl variables are defined in the BLOCK DATA subprogram. 
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/C0NST2/ 


COMMON / CONST2/DPI, DTWOPI, DRAD, DRSEC, PI, TWOPI, RAD, RSEC 
Data base containing constants necessary for angular calculations in radians 


Variable 

Type 

Description 

DPI 

D 

vin radians, 3. 141592653 5897932 

DTWOPI 

D 

2 77 in radians, 6.283 1853 071795864 

DRAD 

D 

Conversion from degrees to radians, .017453 292519943 2 96 

DRSEC 

D 

Conversion from seconds of arc to radians, .484813 681109536 x 10 ' 5 

PI 

R 

77 in radians, 3.141593 

TWOPI 

R 

2t 7 in radians, 6.283185 

RAD 

R 

Conversion from degrees to radians, .01745329 

RSEC 

R 

Conversion from seconds of arc to radians, .4848137 x 10' 5 


The CONST2 variables are defined in the BLOCK DATA subprogram. 
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/ C0NST3/ 


COMMON /CONST3/THDOT2, THDT2S, AESQ, GM3(2), IY1, B, BO, APGM, 

APLM, FSQ32, FFSQ32 


Program 

Where 


Variable 

Type 

Description 

Defined 

THDOT2 

D 

Total rotation of the mean Greenwich meridian 
with respect to the mean equinox of date - radians 
per 24 hours of UTl time 

BLOCK 

DATA 

TKDT2S 

D 

Total rotation of the mean Greenwich meridian 
with respect to the mean equinox of date - radians 
per second of UTl time (same as THDOT2 except 
for units) 

BLOCK 

DATA 

AESQ 

D 

R E 2 - meters 2 

ORBIT 

GM3 (I) 

R 

3GM m - meters 3 /seconds 2 

ORBIT 

GM3 (2) 

R 

3GM 0 - meters 3 /seconds 2 

ORBIT 

IY1 

I 

Years elapsed from 1959 to the year of the epoch 
of the initial elements 

ORBIT 

B 

R 

1^7 C D (= 0 if M s < 0) 

ORBIT 

BO 

R 

1 A s 

2 M s (= 0 if M s < 0) 

ORBIT 

APGM 

R 

p o C R (= 0ifM s < 0) 

ORBIT 

APLM 

R 

IT p r, (= 0 if M s < 0) 

s 

ORBIT 

FSQ32 

R 

yR E f 2 

APPER 

ORBIT 

FFSQ32 

R 

R E f i 1 + 2 f ) 

APPER 


ORBIT 
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/CONST3/ (continued) 


where 

G = universal gravitational constant 

M M = mass of moon in earth masses 

M 0 = mass of sun in earth masses 

M s = mass of satellite - kilograms 

R e = semi-major axis of earth's reference ellipsoid 

f = flattening of earth's reference ellipsoid 

A s = satellite cross sectional area - meters 2 

C D = satellite drag coefficient 

C E = satellite reflectivity coefficient 

P o = solar radiation pressure in the vicinity of the 
earth - newtons/meter 2 

The following variables are defined in BLOCK DATA: 

THDOT2 = 6.300 388 098 445 593 radians per 24 hours of UTl time 
THDT2S = .729 211 585 468 2 X 10" 4 radians per second of UTl time. 

The remaining variables are calculated in subprogram ORBIT. 
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/CONVkG/ 


Variable 
TORE FT 


COMMON / CONVRG/TOREFT 


Type 


Description 


Program 

Where 

Defined 


L Initialized to FALSE in BLOCK BLOCK 

DATA. The OUTPUT option DATA 

card is used to set this switch OPTCRD 

to TRUE which indicates an 
ORB1 tape has been requested 
and also to reference all quanti- 
ties to the reference date. 


Program 

Where 

Used 


MAIN 

ORBIT 


/ C0RB1/ 


Variable 

RANDOT 

PERDOT 

PERHT 

APHT 

PRD 


COMMON /CORB1/RANDOT, PERDOT, PERHT, APHT, PRD 


Type 

Description 

Program 

Where 

Defined 

Prograr 

Where 

Used 

D 

Rate of change of the right 
ascension of ascending node 
(deg/day) 

MAIN 

MAIN 

D 

Rate of change of the argument 
of perigee (deg/day) 

MAIN 

MAIN 

D 

Perigee height, considering 
oblateness (km) 

APPER 

MAIN 

D 

Apogee height, considering 
oblateness (km) 

APPER 

MAIN 

D 

Period (seconds) 

MAIN 

MAIN 
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/ COSAVE/ 


COMMON / COSAVE/SA(5, 9), S2, S3, SCENTE 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

SA(5, 9) 

D 

The saved coefficients of the 
polynomials used to determine 
the ephemeris quantities 

EPHQAN 

FRCS 

S2 

D 

The overlap span from one arc 
of ephemeris quantities to the 
next 

EPHQAN 

EPHQAN 

S3 

D 

The length of arc used to fit the 
ephemeris quantities. The 
nominal value set in the program 
is 2.5 days 

BLOCK 

DATA 

EPHQAN 

SCENTE 

D 

The saved reference time point 
for the ephemeris polynomials 

EPHQAN 

FRCS 
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/cpartl/ 


COMMON /CPARTL/ PXPX0(50, 6), PMPX0(50, 2), PMPSTA0, 2) 


Program 

Where 


Variable 

Type 

Description 

Defined 

PXPXO (50, 6) 

D 

Array containing partials of posi- 

ORBIT 


tion and velocity with respect to 
the parameters to be estimated. 
The matrix is of the form: 


dx 

dy 

dz 

dx 

dy 

d Z 

dx 0 

dx 0 

dx 0 

dx 0 

dx Q 

dx 0 

dx 

dy 

dz 

d x 

dy 

d z 

d y 0 

d y<> 

d Vo 


a y 0 

^y 0 

dx 

dy 

dz 

dx 

dy 

dz 

1 ° 
N 

I*** 

dz 0 

dz 0 

dz Q 

dz 0 

dz Q 

dx 

dy 

d z 

dx 

dy 

d z 

dk 0 

di 0 

di 0 

dk 0 


dk 0 

d.x 

dy 

d z 

dx 

dy 

d z 



*0 

a y 0 

a y 0 

dFo 

dx 

dy 

dz 

dx 

dy 

d z 

dz 0 

dz 0 

d ' z o 

di 0 


dz 0 

dx 

dy 

dz 

dx 

dy 

d z 

5C d 

5 C d 

5c d 

* c d 

5C D 

5c d 

dx 

dy 

d z 

dx 

dy 

dz 

w 

dC R 

Jc R 

dC R 

5C r 


dx 

dy 

d z 

dx 

dy 

d z 

dC 

nm 

dC 

nm 

dC 

nm 

dC 

nm 

dC 

nm 

dC 

nm 


Program 

Where 

Used 

PREDCT 

MAIN 
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/CPARTL/ (continued) 


Variable Type 


Description 


Program 

Where 

Defined 


PMPX0(50, 2) 


D 


The partials of the measurements 
at a given time with respect to 
the parameters being estimated. 
The matrix is of the form: 


PREDCT 
(calculated 
only when M 
and M 2 have 
non-zero 
sigmas) 


<9M X c?M 2 


dx 0 

a x <, 


<?m 2 

d y 0 

d y* 


au 2 

dz Q 




dk 0 

ai 0 


5M 2 

dy 0 


3M 1 

3M 2 

d'z 0 

d'z 0 

aMj 

5M 2 

ac D 

5C d 


5M 2 

5C r 

3C r 


5M 2 


NOTE : The PXPXO and PMPXO arrays are packed so that if a particular parameter 
to be estimated, the lower parameter partials move up. 


Program 

Where 

Used 


MAIN 
(as input 
to subpro 
gram 
ESTIM) 


type is not 
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/ CPARTL/ (continued) 


Variable 
PMPSTA(3 , 2) 


Typ e 


Description 


Program 

Where 

Defined 


D 


The partials of the measurements 
at a given time with respect to 
tracking station coordinates. The 
matrix is of the form: 


<3Mj 

<?m 2 


du l 

<JMj 

3M 2 


<5u 2 


dU 2 


du 3 


PREDCT 
(calculated 
only for 
observations 
with non- 
zero sigmas 
made by 
stations 
whose 
coordinates 
are to be 
adjusted) 


where 



0 


geocentric inertial rectangular 
coordinates of the position and 
velocity vectors of the satellite 
at the initial epoch 


Program 

Where 

Used 

MAIN 
(as input 
to subpro- 
gram 
ESTIM) 
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/CPARTL/ (continued) 


C D = satellite aerodynamic drag coefficient 


C R = satellite reflectivity coefficient 


Mjl _ observation or observations of the satellite which 
M 2 j is one of the following types or pair of types: 



right ascension 
declination 


[Mj range 
< Mj >= range rate 

rM J range rate (derived from Doppler) 



minitrack direction cosine 
minitrack m direction cosine 



X angle (30 foot antenna) 
Y angle (30 foot antenna) 



aaimuth angle 
elevation angle 


uA geocentric earth-fixed rectangular 
u 2 j= coordinates of the tracking station 
u 3 J (computed from its geodetic spherical 

coordinates); u 3 axis is directed 
toward earth's positive direction of 
axis of rotation; u x axis is directed 
toward intersection of Greenwich 
meridian and earth's equator. 
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/CQUANT/ 


Variable 

STAX(N) 

STAT(N) 

STAZ(N) 


NHAT(N) 

ZHAT(N) 

EHAT(N) 


THPRIM(l, N) 


THPRIM(2, N) 


COMMON / CQUANT/ ST AX (50), STAY(50), STAZ(50), NHAT(3, 50), 
ZHAT(3, 50), EHAT(3, 50), THPRIM(2, 50) 



D 


sin (<£ N ~<£ N ') 


D 


cos (4 > n - 0 n *) 


Program 

Program 

Where 

Where 

Defined 

Used 

SQUANT 

MAIN 

(Calcu- 

ESTIM 

lated on 

OBSDOT 

1st call 

PLHOUT 

only) 

PREDCT 

SQUANT 

(Calcu- 
lated on 

OBSDOT 

1st call 
only) 

PREDCT 

SQUANT 

Not cur- 
rently 
used by 
GEOSTAR 
I 

SQUANT 

Not cur- 
rently 
used by 
GEOSTAR 


I 


where: 


4> u = geodetic latitude of Nth station 
tp N ' = geocentric latitude of Nth station 



geocentric earth fixed rectangular coordinates of the Nth 
tracking station (computed from its geodetic spherical co- 
ordinates); u 3 axis is directed towards earth's positive 
direction of axis of rotation; axis is directed towards 
intersection of Greenwich meridian and earth’s equator 
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/ C QUANT/ (continued) 


unit vectors centered at the Nth tracking 
_ stations directed North, Vertical and East 
respectively. The components of these 
vectors are referred to the topocentric 
equatorial, earth fixed coordinate system 
(principal axis parallel to the principal 
axis of the geocentric system). 



/CSTHET/ 


COMMON / CSTHET/ CTHETG, STHETG 





Program 

Where 

Variable 

Type 

Description 

Defined 

CTHETG 

D 

Cosine of the apparent right 
ascension of the mean Greenwich 
meridian 

PREDCT 

STHETG 

D 

Sine of the apparent right as- 
cension of the mean Greenwich 
meridian 

PREDCT 


Program 

Where 

Used 

XEFIX 


XEFJX 
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/CSTENF/ 

COMMON / CSTINF/MEASNQ{4), NOBS(4), RDMEAN(4), RMS0(4), RND(4), 
MEASWT(4), WTMEAN(4), RMSWTQ(4), WTRND(4), TYPRMS(14), NOTYPE(14) 

NOTE 

For each data reduction and parameter estimation 
run made, the following restrictions apply: 

50 = maximum number of 

tracking stations allowed 

4 = maximum number of dif- 
ferent observation types 
allowed per station.. 

At the end of each iteration, subprogram STAINF is 
called to provide a statistical summary for that 
particular iteration on a station by station basis. 

At the end of the iteration, STAINF calculates the 
required statistical information for each of the 
observation types (maximum 4) for each station 
and returns this information to MAIN via 
COMMON/ CST INF/'. MAIN then prints this 
statistical summary for each station. The di- 
mension of 4 shown in the above COMMON/ 

CSTINF/ block for each variable represents the 
maximum number of observation types allowed 
per station. At any given time, the information 
in this block pertains to only one station. The statis- 
tical information corresponding to the same in- 
dex in the various arrays pertains to the specific 
observation type implied by that particular index 
number. 
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/CSTESTF/ (continued) 


Variable Type Description 


Program 

Where 

Defined 


NOTE 


Let N (= 1, 2, 3, or 4) be the index in the array 
indicating a specific observation under considera- 
tion as described in the following table : 

MEASNO(N) I Code number (MTYPE) indicating STAINF 

type of observation for which sta- 
tistical information was computed. 



MTYPE 

= 1 

for right ascension 




= 2 

for range 




= 3 

for range rate 




= 4 

for range rate (derived 





from Doppler) 




= 5 

for direction cosine 




= 6 

for X angle 




= 7 

for azimuth angle 




= 8 

for declination 




= 12 

for m direction cosine 




= 13 

for Y angle 




= 14 

for elevation angle 




= 16 

for geocentric coordinate x 




= 17 

for geocentric coordinate y 




= 18 

for geocentric coordinate z 




» 19 

for geocentric coordinate x 




= 20 

for geocentric coordinate y 




= 21 

for geocentric coordinate z 


NOBS(N) 

I 

Total number of MTYPE observa- 

STAINF 



tions for this station 


RDMEAN(N) 

R 

Mean of the residuals of all 

STAINF' 



MTYPE observations for this 




station 



Program 

Where 

Used 

MAIN 


MAIN 

MAIN 
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/CSTINF/ (continued) 


Variable 

Type 

Description 

Program 

Where 

Defined 

Pi ogram 
Where 
Used 

RMSO(N) 

R 

RMS of the residuals about zero 
for all MTYPE observations for 
this station (calculated only if 
there are 10 or more MTYPE 
observations for this station) 

STAINF 

MAIN 

RND(N) 

R 

Random normal deviate of resi- 
duals for all MTYPE observa- 
tions for this station (calculated 
only if there are 10 or more 
MTYPE observations for this 
station) 

STAINF 

MAIN 

MEASWT(N) 

I 

Number of weighted MTYPE ob- 
servations for this station 

STAINF 

MAIN 

WTMEAN(N) 

R 

Weighted mean of the residuals 
for all weighted MTYPE obser- 
vations for this station 

STAINF 


RMSWT0(N) 

R 

Weighted RMS of the residuals 
about 0 for the weighted MTYPE 
observations for this station 
(calculated only if there are 10 
or more observations for this 
station) 

STAINF 

MAIN 

WTRND(N) 

R 

Weighted random normal deviate 
of the residuals of the MTYPE 
observations for this station 
(calculated only if there are 10 
or more observations for this 
station) 

STAINF 

MAIN 

TYPRMS 

R 

Total RMS for each observation 
type in the solution 

STAINF 

MAIN 

NOTYPE 

I 

Number of observations per ob- 
servation type if zero, that data 
type is not in A he solution 

STAINF 

MAIN 
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/ csvevl/ 


COMMON / CSVEVL/ CNM(3, 3), SNM(3, 3) 

These arrays are used only by subroutine VEVAL. They are initialized by STORGP from the data 
base array of geopotential coefficients /FMODEL/CS. 


Variable Type 


Description 


CNM(3, 3) 
SNM(3 , 3) 


D Contain constants for the geopotential coefficients C or S; zonals 

D through 4th order and tesserals through 3rd order are included. 


Variable Name 
with FORTRAN Subscript 

CNM(1, 1) 

CNM(2, 1) 

CNM(3, 1) 

CNM(1, 2) 

CNM(2, 2) 

CNM(3, 2) 

CNM(1, 3) 

CNM(2, 3) 

CNM(3, 3) 

Same as above for SNM 


Values Using 
Geophysical Notation C n 

C(2, 0) 

C(2, 1) 

C(3, 1) 

C(3, 0) 

C(2, 2) 

C(3, 2) 

C(4, 0) 

C(3, 3) 
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/CVRCOV/ 


COMMON / CVRCOV/VRCOV (50, 50), JHIGH, VRCSW 





Program 

Program 




Where 

Where 

Variable 

Type 

Description 

Defined 

Used 

VRCOV 

R 

Initial variance -covariance 

BLOCK 

ESTIM 



matrix, the inverse of which 

DATA or 

MAIN 



is the parameter weight matrix 

input on 




used in ESTIM. BLOCK DATA 

cards in 




defines the six diagonal ele- 

OPTCRD 



ments to be used for the state 
vector. The entire matrix may 
be filled by using the VARCOV 
option card. 


JHIGH 


VRCSW 


Initialized as 6. If VARCOV 
option cards are read in, JHIGH 
indicates the dimension of the 
portion of VRCOV that is filled. 

BLOCK 

DATA 

OPTCRD 

ESTIM 

MAIN 

Initialized to FALSE; if VARCOV 
option cards are read in, VRCSW 
is automatically set to TRUE in- 
dicating that the initial variance- 

BLOCK 

DATA 

OPTCRD 

ESTIM 


covariance matrix has been 
redefined. 
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/DC/ 

COMMON /DC/H(6, 6), DETER, MMl 


Variable Type 


Description 


Program 

Where 

Defined 


H(6, 6) 


DETER 

MMl 


D Matrix containing the Jacobian CSTEP 

of accelerations with respect 
to position and velocity 

D Determinant of matrix inverted INV2 

by INV2 

I Dimension of matrix to be CSTEP 

inverted 


Program 

Where 

Used 


INV2 


INV2 
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/DRGBliK/ 


COMMON /DRGBLK/ C (4) , SPSISQ, C3, Cl, VEL, XDOTR, YDOTR, RHO, HT 





Program 

Where 

Program 

Where 

Variable Type 

Description 

Defined 

Used 

C(l) " 





C(2) 

C(3) 

C(4) 

J 

>• R 

Used for transfer of intermediate 
calculations 

DENSTY 

VEVAL 

SPSISQ 

R 

z 2 

— r = sin 2 4> 

DRAG 

VEVAL 

VEVAL 

C3 

BBRHO . 

\ 

P 

V r 

DRAG 

VEVAL 

VEVAL 

Cl 

BBRHO V. 

1 

1 p A s 

2 V r M s - a D 

DRAG 

VEVAL 

VEVAL 

VEL 

R 

v r - meters/second 

DRAG 

VEVAL 

VEVAL 

XDOTR ' 
VELR(l) . 

\ 

x r - meters/second 

DRAG 

VEVAL 

VEVAL 

YDOTR ] 
VELR(2) , 

\ 

y r - meters/second 

DRAG 

VEVAL 

VEVAL 

RHO 

VELR(3) . 

\ 

p - kilograms/meter 3 

DRAG 

VEVAL 

VEVAL 

HT 

R 

h - meters 

DRAG 

VEVAL 

DENSTY 

VEVAL 


where 

z = inertial z component of satellite 
position vector 

r = geocentric distance of satellite 

p = atmospheric density at the height 
of the satellite 
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/DRGBLK/ (continued) 


A s = satellite cross-sectional area- meters 2 
M s = satellite mass - kilograms 
C D = satellite drag coefficient 
v r = magnitude of relative velocity 
h = height of satellite 

x r = inertial x component of relative velocity 
vector of satellite 

y r = inertial y component of relative velocity 
vector of satellite 

<P = geocentric latitude of satellite 

a D = magnitude of acceleration due to aerodynamic 
drag 

The variables C(l), C(2), C(3), and C(4) are coefficients of the equation 


log l0 P = C(l) + C(2)h + C(3)h 2 + C(4)h 3 . 

£ 

Subroutine VEVAL requires the quantity^ iog 10 p , so that only the coefficients C(2), C(3) and 
C(4) are used by that program. The corresponding variables in VEVAL are: 


DENSTY 

VEVAL 

C(l) 

CO 

C(2) 

EXPT(l) 

C(3) 

EXPT (2) 

C(4) 

EXPT (3) 

XDOTR = 

VELR(l) 

YDOTR = 

VELR(2) 

RHO = 

VELR(3). 


After calculation using these variables, VEVAL then uses CO and EXPT as scratch storage. The 
equivalence of the same variables in DRAG and VEVAL with different names are shown in the above 
list in brackets. These are 
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C3 = BBRHO 


Cl = BBRHOV 



/ESTGP/ 


COMMON /ESTGP/ CSA(50), CSE(50), CSSIG(50), PERTRB, LABELG, NCSN(50), 
NCSM(50), NC, NS, NPRTL, ITERGP, 1ST ATE 

/ESTGP/ defines variables used when estimating the earth's geopotential coefficients or gravity 
parameters (GP). 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

CSA(50) 

R 

Array of a priori values of 
geopotential coefficients for 
estimation 

At first 
iteration 
CSE = CSA. 

ESTIM 

CSE(50) 

R 

Array of estimated values of 
geopotential coefficients for 
estimation 

READGP 
from Data 
Base, 

/FMODEL/ 
CS(30, 33) or 
input cards. 

ESTIM 
MAIN . 
WTBMAT 

CSSIG(50) 

R 

Array of a priori standard 
deviations of the geopotential 
coefficients for estimation 

READGP 

ESTIM 

PERTRB 

R 

Not used 



LABELG 

I 

Position in the GRPAR array 
corresponding to the gravity 
coefficients. LABELG is 
initialized to zero in BLOCK 
DATA. 

BLOCK 

DATA 

OPTCRD 

EGRAV 
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/ESTGP/ (continued) 


Program Program 

Where Where 


Variable 

Type 

Description 

Defined 

Used 

NCSN(50) \ 

I 

Arrays of n and m subscripts 

Input cards 

MAIN 

NCSM(50) J 


designating which C nm or S,, m 

as proc- 

EGRAV 



are represented in CSA and 

cessed by 

WTBMAT 



CSE. NCSN contains the n 

READGP 

STATRD 


subscript; NCSM, the m sub- 
script. Subscripts for C nm 
are first in the arrays, with 
S nm following. The first NC 
places are for C nm and the 
following NS places are for 


n,m • 


NC 

I 

Number of C nm to estimate 

READGP 

MAIN 

EGRAV 

STATRD 

WTBMAT 

NS 

I 

Number of S nm to estimate 

READGP 

MAIN 

EGRAV 

STATRD 

WTBMAT 

NPRTL 

I 

Total number of GP to esti- 
mate, NC + NS which is 
initialized to 0 in BLOCK 
DATA 

BLOCK 

DATA 

READGP 

MAIN 

OPTCRD 

ORBIT 

STATRD 

ESTIM 

EGRAV 

VEVAL 

WTBMAT 

ITERGP 

I 

Number of iterations to com- 
pute numerical partials; de- 
fined as 0 in BLOCK DATA and 

BLOCK 

DATA 

READGP 

MAIN 


redefined on the COEFGP option 
card 
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/ESTGP/ (continued) 


Variable 

ISTATE 


Type 


Description 


Program 

Where 

Defined 


Program 

Where 

Used 


L 


The switch is initialized to BLOCK 

TRUE in BLOCK DATA, indi- DATA 

eating that state parameters READGP 

are to be estimated. The 
COEFGP option card is used to 
set this switch to FALSE which 
indicates that state parameters 
will not be estimated. 


MAIN 

STATRD 

ORBIT 

SWTEST 

WTBMAT 
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/FLXBLK/ 


COMMON / FLXBLK/DSTART, DAY, AVFLX(15), DFLX(15), AP(15) 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

DSTART 

D 

Epoch of initial elements - days 
from January 0.0 of the year of 
the epoch of the initial elements 

MAIN 

DENSTY 

GEOSRD 

PREDCT 

DAY2 

D 

Epoch of calculation - days from 
January 0.0 of the year of the 
epoch of the initial elements 

FRCS 

DENSTY 

AVFLX(15) 

R 

Solar flux values averaged over 
55 days: 

MAIN 

DENSTY 

DFLX(15) 

R 

Daily values of the solar flux for 
15 days starting with DSTART 

MAIN 

DENSTY 

AP(15) 

R 

Geomagnetic activity index 
(dimensionless) 

MAIN 

DENSTY 


Units of AVFLX(15) and DFLX(15) are 1<T 22 
watts/meter 2 /cycle/sec of bandwidth. 
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/FMODEL/ 


COMMON / FMODE L/ MODE L (8) , CS(30, 33), INDEXl, INDEX3 


Program Program 

Where Where 


Variable 

Type 

Description 

Defined 

Used 

MODEL (8) 

D 

Alphanumeric information describ- 
ing gravity coefficient set used. 
This array is printed by MAIN 
for identification purposes. 

BLOCK 

DATA 

MAIN 

COEFL 

EGRAV 

CS(30, 33) 

R 

The C and S coefficients in the 
spherical harmonic expansion 
of the expression for the 
geopotential 

BLOCK 

DATA 

MAM 

COEFL 

EGRAV 

READGP 

INDEXl 

I 

N = index of highest degree 
of spherical harmonic 
coefficients contained in 
CS 

BLOCK 

DATA 

MAIN 

COEFL 

EGRAV 

INDEX3 

I 

M = index of highest order of 
spherical harmonic coef- 
ficient contained in 

BLOCK 

DATA 

MAIN 

COEFL 

EGRAV 


CS 

The C and S coefficients are described by indices 
n and m. For the zonal harmonics, m = 0. For 
the sectorial harmonics, m = n. For the tesseral 
harmonics, m < n. Therefore, the C coefficients 
only fill one triangle of a matrix, as do the S. To 
conserve space, both sets of coefficients were com- 
bined into one matrix. A diagram of that matrix 
and the computations for the subscripts correspond- 
ing to the n and m indices appear in the EGRAV 
subroutine description in Section 4.2. The maximum 
geopotential field allowed is (30x30). A(15xl5) geo- 
potential (modified SAO C - 5 gravity model) is 
currently defined in BLOCK DATA. The coefficients 
of this model are listed in Table 3. 
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/FMODEL/ (continued) 


C(2, 0) 
S(2, 0) 
C (2, 2) 
S(2, 2) 
C(3, 0) 
S(3, 0) 
C(3, 1) 
S(3, 1) 
C(3, 2) 
S(3, 2) 
C(3, 3) 
S(3, 3) 
C(4, 0) 
S(4, 0) 
C(4, 1) 
8(4, 1) 
C(4, 2) 
S(4, 2) 
C(4, 3) 
S(4, 3) 
C(4, 4) 
S(4, 4) 
C(5, 0) 
S(5, 0) 
C( 5, 1) 
S(5. 1) 
C(5, 2) 
S(5, 2) 
C(5, 3) 
S(5, 3) 
C(5, 4) 
S(5, 4) 
C(5, 5) 
S(5, 5) 
C(6, 0) 
S(6, 0) 
C(6, 1) 
S(6, 1) 
C(6, 2) 
S(6, 2) 
C{6, 3) 
S(6, 3) 
C(6, 4) 
S(6, 4) 
C(6, 5) 
S(6, 5) 


Table 3 

SAO Denormalized Coefficients (C -• 5) 


= - 1082.645 • 10" 6 

C ( 6 , 6 ) 

= 

- 0 . 932 - 

10* 11 

C ( 10 , 02 ) 


- 0.624 • 

10 " 8 

= 0.0 


S ( 6 , 6 ) 

= 

- 0.361 • 

10" 10 

S ( 10 , 02 ) 

= 

- 0.250 • 

10 " 8 

= + 1.536 • 

10 " 6 

C ( 7 , 0 ) 

= 

+ 0.333 • 

10" 6 

C ( 10 , 03 ) 


- 0.379 • 

10 " 9 

= - 0 . 872 - 

10" 6 

S ( 7 , 0 ) 

= 

0.0 


S ( 10 , 03 ) 


+ 0.175 • 

10 " 9 

= + 2 . 546 - 

10' 6 

C ( 7 , 1 ) 


+ 0 . 144 - 

10" 6 

C ( 10 , 04 ) 

= 

- 0.436 • 

10' 10 

= 0.0 


S ( 7 , 1 ) 

= 

+ 0.114 • 

10" 6 

S ( 10 , 04 ) 


- 0.654 • 

10' 10 

= + 2 . 091 - 

10' 6 

C ( 7 , 2 ) 

= 

+ 0.363 • 

10" 7 

C ( ll , 0 ) 

= 

- 0.302 • 

10~ 6 

= + 0.287 • 

10' 6 

S ( 7 , 2 ) 

= 

+ 0.162 • 

10" 7 

8 ( 11 , 0 ) 

= 

0.0 


= + 0.251 • 

10' 6 

C ( 7 , 3 ) 

= 

+ 0.352 • 

10" 8 

C ( ll , 01 ) 

= 

- 0.313 • 

10- 7 

= - 0 . 184 - 

10' 6 

S ( 7 , 3 ) 

= 

+ 0.254 ■ 

io' 9 

S ( ll , 01 ) 

=r 

+ 0.880 • 

10' 8 

= + 0.782 • 

10~ 7 

C ( 7 , 4 ) 

= 

- 0.323 • 

10" 9 

C ( 12 , 0 ) 

= 

+ 0.357 • 

10" 6 

= + 0.226 • 

10" 6 

S ( 7 , 4 ) 

=r 

- 0.217 • 

10" 9 

S ( 12 , 0 ) 

= 

0.0 


= + 1.649 • 

10 “ 8 

C ( 7 , 5 ) 

= 

+ 0.269 ■ 

IO" 10 

C ( 12 , 01 ) 

= 

- 0.400 • 

10' 7 

= 0.0 


S ( 7 , 5 ) 

= 

+ 0.191 • 

10' 10 

S ( 12 , 01 ) 

= 

- 0.402 • 

10" 7 

= - 0.543 • 

10' 6 

C ( 7 , 6 ) 

= 

- 0.145 • 

10- 10 

C ( 12 , 02 ) 

= 

- 0.470 • 

10' 8 

= - 0.445 • 

10~ 6 

S ( 7 , 6 ) 

= 

+ 0.437 • 

10 -11 

S ( 12 , 02 ) 

=: 

- 0.230 » 

10" 9 

= + 0.738 • 

10" 7 

C ( 7 , 7 ) 

= 

+ 0.102 • 

10 -11 

C ( 12 , 12 ) 

- 

- 0.278 • 

10" 18 

= + 0.148 ■ 

10" 6 

S ( 7 , 7 ) 

= 

+ 0.180 • 

10" 11 

S ( 12 , 12 ) 

= 

+ 0.718 ■ 

10-20 

= + 0.509 • 

10 -7 

C ( 8 , 0 ) 

= 

+ 0.270 • 

10" 6 

C ( 13 , 0 ) 

= 

+ 0.114 ■ 

10' 6 

= - 0.114 • 

10* 7 

S ( 8 , 0 ) 

= 

0.0 


S ( 13 , 0 ) 

= 

0.0 


= - 0.112 • 

10~ 8 

C ( 8 , 1 ) 

= 

- 0.520 • 

10" 7 

C ( 13 , 12 ) 

= 

- 0.126 • 

10" 18 

= + 0.486 • 

10' 8 

8 ( 8 , 1 ) 

= 

+ 0.450 • 

10" 7 

S ( 13 , 12 ) 

= 

+ 0.117 • 

10" 18 

= + 0.210 • 

10" 6 

C ( 8 , 2 ) 

= 

+ 0.214 ■ 

10" 8 

C ( 13 , 13 ) 

= 

- 0.239 • 

10 -19 

= 0.0 


8 ( 8 , 2 ) 

= 

+ 0.320 • 

10 -8 

S ( 13 , 13 ) 

- 

+ 0.212 • 

10‘ 19 

= - 0.677 • 

10- 7 

C ( 8 , 3 ) 

= 

- 0.374 • 

10' 9 

C ( 14 , 0 ) 

= 

- 0.179 • 

10" 6 

= - 0.882 • 

10' 7 

S ( 8 , 3 ) 

= 

+ 0.404 • 

10’ 10 

8 ( 14 , 0 ) 

= 

0.0 


= + 0.102 • 

10 -6 

C ( 8 , 4 ) 

= 

- 0.277 • 

10~ 9 

C ( 14 , 01 ) 

= 

- 0.788 • 

10 " 8 

= - 0.375 • 

10" 7 

8 ( 8 , 4 ) 

= 

- 0.157 • 

10' 10 

S ( 14 , 01 ) 

= 

+ 0.280 ■ 

10 " 8 

= - 0.172 • 

10" 7 

C ( 8 , 5 ) 

= 

- 0.959 • 

10" 11 

C ( 14 , 11 ) 

ss 

+ 0.947 • 

10 " 21 

= + 0.231 • 

10~ 9 

8 ( 8 , 5 ) 

= 

+ 0.214 • 

10 _1 ° 

S ( 14 , 11 ) 

= 

- 0.473 • 

10- 21 

= - 0.206 • 

10" 8 

C ( 8 , 6 ) 

= 

- 0.475 • 

10- 12 

C ( 14 , 12 ) 

= 

+ 0.140 - 

10 ” 2 0 

= + 0.498 • 

10" 9 

8 ( 8 , 6 ) 


+ 0.888 • 

10- 11 

S ( 14 , 12 ) 

= 

- 0.132 • 

10" 19 

= + 0.384 • 

10" 9 

C ( 8 , 7 ) 

= 

- 0.444 • 

10 -13 

C ( 14 , 14 ) 

= 

- 0.193 • 

10' 21 

= - 0.146 • 

10' 8 

8 ( 8 , 7 ) 

= 

+ 0.158 ■ 

10" 12 

S ( 14 , 14 ) 

= 

- 0.414 • 

10' 22 

= - 0.646 • 

10" 6 

C ( 8 , 8 ) 

= 

- 0.316 • 

10 -12 

C ( 15 , 09 ) 

= 

- 0.241 • 

10- 18 

= 0.0 


8 ( 8 , 8 ) 

= 

+ 0.130 • 

10” 12 

S ( 15 , 09 ) 

= 

- 0.483 • 

10' 18 

= - 0.370 ■ 

10- 7 

C ( 9 , 0 ) 

= 

+ 0.530 • 

10' 7 

C ( 15 , 12 ) 

=. 

- 0.138 • 

10" 19 

= - 0.212 • 

10' 7 

S ( 9 , 0 ) 

= 

0.0 


S ( 15 , 12 ) 


- 0.190 • 

10 " 2 ° 

= + 0.858 • 

10' 8 

C ( 9 , 1 ) 

= 

+ 0.760 • 

10" 7 

C ( 15 s 13 ) 

= 

- 0.770 • 

10~ 21 

= - 0.455 • 

10" 7 

S ( 9 , 1 ) 

= 

+ 0.784 • 

10 " 8 

S ( 15 , 13 ) 

= 

- 0.374 • 

10" 21 

= - 0.112 • 

10" 8 

C ( 9 , 2 ) 

= 

- 0.277 • 

10' 9 

C ( 15 , 14 ), 

= 

+ 0.114 • 

10- 2 2 

= + 0.643 • 

10" 9 

8 ( 9 , 2 ) 

= 

+ 0.242 • 

10" 8 

S ( 15 , 14 ) 

= 

- 0.558 • 

10~ 22 

- - 0.167 • 

10- 9 

C ( 10 , 0 ) 

= 

+ 0.540 • 

10 -7 

C ( 17 , 13 ) 

= 

+ 0.159 • 

10- 22 

= - 0.196 • 

10~ 8 

S ( 10 , 0 ) 

= 

0.0 


S ( 17 , 13 ) 


+ 0.280 ■ 

10" 22 

= - 0.253 • 

10' 9 

C ( 10 , 01 ) 

= 

+ 0.649 • 

10 -7 





= - 0.370 • 

10- 9 

S ( 10 , 01 ) 

= 

- 0.779 - 

10 -7 






126 





Variable 
GRP Alt (3, 50) 


/GRBLOK/ 

COMMON /GRBLOK/ GRPAR(3, 50) 


Type 


Description 


Program Program 

Where Where 

Defined Used 


D 


Array of acceleration partials 
with respect to the parameters 
being estimated. The array is 
of the form: 


<9x 

dx 

dx 

3C d 

ac R 

ac 

nm 

dy 

dy 

dy 

dC B 

dC R 

ac 

nm 

dz 

dz 

dz 

dC D 

dC R 

dC 

nm 


DRAG 

FRCS 

EGRAV 

VEVAL 


ORBIT 

RK 

CSTEP 

TABLEB 


where 





coordinates of the satellite acceleration 


and 


C D = drag coefficient, 
c R = reflectivity coefficient, 
c nm = geopotential coefficient. 
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/IORBIT/ 


COMMON / IORBIT/ TREQ, NTER, INITSW 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

TREQ 

D 

Observation time (or output 
time) 

MAIN 

ORBIT 

HEMINT 

NTER 

I 

Iteration number which is 
initialized to 1 in BLOCK 
DATA 

BLOCK 

DATA 

MAIN 

ORBIT 

INITSW 

L 

Switch to initialize the orbit 
generator 

MAIN 

ORBIT 

EGRAV 

ORBl 
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/LIMITS/ 

COMMON /LIMITS/S 1(3, 50), S2(3, 50), CTOL, ITER, ISCT(2), 
ISWT(30), SW(30), ORDER (2) 


Program Program 





Where 

Where 

Variable 

Type 

Description 

Defined 

Used 

SI (3, 50) 

D 

First sums necessary for the 

SUMS 

CSTEP 



predictor-corrector formulas 


TABLEB 

S2(3, 50) 

D 

Associated second sums 

SUMS 

CSTEP 

TABLEB 

CTOL 

D 

Predictor- corrector tolerance 

BLOCK 

DATA 

OPTCRD 

CSTEP 

ITER 

I 

Number of predictor-corrector 
iterations required to satisfy 
CTOL tolerance (must be <3) 

CSTEP 

ORBIT 

ISCT(2) 

I 

ISCT(l)-number of points pro- 

ORBIT 

TEST 



duced at a particular stepsize 
for equations of motion 

ISCT(2)-number of points pro- 
duced at a particular stepsize 
for the variational equations 


TABLEB 

ISWT (3 0)* 

L 

Array of logical external switches 

BLOCK 

FRCS 



which set various options in the 

DATA 

ORBIT 



program 

OPTCRD 

SWTEST 




MAIN 

STATRD 

SW(30)* 

L 

Array of logical internal switches 

BLOCK 

EGRAV 



which are set on condition during 

DATA 

VEVAL 



the execution of a run 

ORBIT 

CSTEP 




FRCS 

HEMINT 


FRCS 

SWTEST 

TABLE 

TABLEB 

TEST 
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/LIMITS/ (continued) 


Variable 


ORDEIt(2) 




Program 

Program 



Where 

Where 

Type 

Description 

Defined 

Used 

I 

ORDER(l)-order of the multistep 

BLOCK 

HEMINT 


predictor -corrector integration 

DATA 

TABLE B 


formulas for the equations of 
motion 

ORBIT 

TEST 


ORDER(2)-order of the muitistep 
corrector integration formulas 
for the variational equations 




The following values, as initially assigned by BLOCK DATA, may be changed by optional input 
cards: 


CTOL = 5.0 X 10' 5 
ORDER (1) = 11 
ORDER(2) = 7 


NOTE: Only ISWT numbers (5), (6), (7), (8), (9), (12), (13), (15), (21), (22), (23), (24), (25), (28), (29), 
(30), and SW numbers (1), (2), (3), (6), (8), (16), (20), (22), and (24) are being used. 

Those ISWT and SW numbers which are presently being used are defined in Tables 4 and 5 
respectively. 
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Table '4 


Variable 

Description 

ISWT<5) 

Setting this switch to TRUE means use two-body gravity model 
for variational equations computation. 

ISWT(6) 

Setting this switch to TRUE means omit drag model for 
variational equations computation. 

ISWT(7) 

Setting this switch to TRUE means omit solar radiation 
model for variational equations computation. 

ISWT(8) 

Setting this switch to TRUE means omit solar gravity model 
for variational equations computation. 

ISWT(9) 

Setting this switch to TRUE means omit lunar gravity model 
for variational equations computation. 

ISWT (12) 

Setting this switch to TRUE means use only two-body gravity 
model in position acceleration computation. 

ISWT (13) 

Setting this switch to TRUE means output position partials 
at each data point. 

ISWT (15) 

Setting this switch to TRUE means output Keplerian elements 
at each data point. 

ISWT (21) 

Setting the switch to TRUE means use two-body gravity model 
in position acceleration computation. 

ISWT (22) 

Setting the switch to TRUE means use lunar gravity model 
in position acceleration computation. 

ISWT(23) 

Setting the switch to TRUE means use solar gravity model 
is position acceleration computation. 

ISWT (24) 

Setting the switch to TRUE means use drag model in 
position acceleration computations. 

ISWT (25) 

Setting the switch to TRUE means use solar radiation 
model in position acceleration computations. 

ISWT (2 8) 

Setting the switch to TRUE means a B matrix (normal 
equations) to be created during first iteration. 

ISWT (2 9) 

Setting the switch to TRUE means read position partials 
from disk. 

ISWT(30) 

Setting the switch to TRUE means write position partials 
onto disk. 
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Variable 

SW(.l) 

SW(2) 

SW(3) 

SW(6) 
SW(8) 
SW(16) 
SW (20) 

SW (22) 
SW (24) 


Table 5 


Description 


Setting this switch to TRUE means the initial table values have 
been completed for the equations of motion. 

Setting this switch to TRUE means the initial table values have 
been completed for the variational equations. 

Setting this switch to TRUE means the local error of the table 
values is being tested for the possibility of doubling or halving 
the initial stepsize. 

Setting this switch to TRUE means doubling of stepsize in TABLE 
or a step increase in CSTEP will take place. 

Setting this switch to TRUE means halving stepsize in TABLE or 
a step decrease in CSTEP will take place. 

Setting this switch to TRUE means reset saved ephemeris 
quantities to original values. 

Setting this switch to TRUE means optimize the stepsize at the first 
CSTEP point for the aquations of motion and the variational 
equations. 

Setting this switch to TRUE means integrate the equations of 
motion and the variational equations with the same stepsize. 

Setting this switch to TRUE means either the equations of 
motion or the variational equations have been integrated 
up to the data point— go back and integrate the other to that 
point. 
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/MOONGR/ 



COMMON /MOONGR/RHOM(3 , 2), RHOSQ(2), RH03(2) 

Program 

Where 

Program 

Where 

Variable 

Type 

Description 

Defined 

Used 

RHOM(3, 2) 

D 

RHOM(I, 1) = x - x m 
RHOM(I, 2) = x - x s 

SUNGRV 

VEVAL 

VEVAL 

RHOSQ{2) 

D 

Pi 2 = R m 2 + R2 - 2(x • x m ) 
p 2 2 = R 2 + R 2 - 2(x • x s ) 

SUNGRV 

VEVAL 

VEVAL 

RH03 (2) 

D 

P/ 

P/ 

SUNGRV 

VEVAL 

VEVAL 


where 


x = satellite position vector 

x ra = lunar position vector 

x s = solar position vector 

R = satellite position vector magnitude 

S m = lunar position vector magnitude 

R s = solar position vector magnitude. 
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/N0N1/ 


COMMON /NON1/ORBTSW, GRAVSW, XYZFSW, XYZLSW, PLTLSW, PLHSW, 
TOREFO, STAPSW, DODSTP, PTAPE, ITRUOR, IBTAPE, JTEMP, DATP 


Variable 


ORBTSW 


GRAVSW 


XYZFSW 


XYZLSW 


PLTLSW 


PHLSW 


Type 


Description 


Program 

Where 

Defined 


L 


L 


L 


L 


L 


L 


Setting this switch to TRUE 
means an orbit generator run 
will be made. 

Setting this switch to TRUE 
means that the gravity field 
being used will be printed. 

Setting this switch to TRUE 
means that the ground track 
will be printed only on the first 
iteration of a data reduction 
run. 

Setting this switch to TRUE 
means that the ground track 
will be printed only on the last 
iteration of a data reduction 
rim. 

Setting this switch to TRUE 
means that a GEORGE tape 
has been requested at the end 
of the run. 

A switch which is initialized to 
FALSE indicating that a 
station's position is in rec- 
tangular coordinates and is to be 
stored in STASIG array. Setting 
this switch to TRUE means that 
a station's position is in geodetic 
coordinates and to be stored in 
PLHSIG array. 


BLOCK 
DATA 
OPT CRD 

BLOCK 
DATA 
OPT CRD 

BLOCK 

DATA 

OPTCRD 


BLOCK 

DATA 

OPTCRD 


BLOCK 

DATA 

OPTCRD 


BLOCK 

DATA 

OPTCRD 


Program 

Where 

Used 


MAIN 

STATRD 

OUTPUT 

MAIN 

STATRD 


MAIN 

STATRD 


MAIN 

STATRD 


MAIN 

STATRD 


MAIN 

STATRD 
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/N0N1/ (continued) 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

TOREFO 

L 

Setting this switch to TRUE 
means that an ORBl TAPE 
will be created. 

BX/OCK 

DATA 

OPTCRD 

MAIN 

STATRD 

STAPSW 

L 

A switch which is initialized to 
TRUE. Setting this switch to 
FALSE means that a STAPOS 
option card and station position 
cards are to be read. 

BLOCK 

DATA 

OPTCRD 

MAIN 

STATRD 

DODSTP 

L 

Setting this switch to TRUE 
means that a DOD's formatted 
observation tape is to be used 
during a data reduction run. 

BLOCK 

DATA 

MAIN 

MAIN 

STATRD 

PTAPE 

L 

Setting this switch to TRUE 
means that a simulated data 
tape comprised of rectangular 
coordinate type data will be 
created during an orbit 
generator run. 

BLOCK 

DATA 

OPTCRD 

main 

STATRD 

ITRUOR 

I 

A variable which defines a 
device to input or output 
simulated observational type 
data. This variable is initialized 
to 36. 

BLOCK 

DATA 

OPTCRD 


IBTAPE 

T 

A variable which defines a device 
to output a B matrix (normal 
equations). This variable is 

BLOCK 

DATA 

OPTCRD 

MAIN 


initialized to 18 . 
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/NON!/ (continued) 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

JTEMP 

I 

A variable which defines the type 
of data reduction run that is being 
made. The variable is initialized 

BLOCK 

DATA 

OPTCRD 

MAIN 

STATED 


to 0 which means that a data re- 
duction run using real data will be 
made. 

JTEMP can be redefined to have 
a value of 1 which indicates that a 
data reduction run using simulated 
data will be made. 

JTEMP can also be redefined to 
have a value of 2 which indicates 
that a data reduction run using real 
data will be made and that the 
computed observations will be 
stored on tape, creating a simu- 
lated data tape. 

DATP I A variable which defines a device BLOCK MAIN 

to store and then use real or DATA 

simulated data. OPTCRD 

The variable is initialized to 35 
indicating that real data will be 
stored and used in a subsequent 
iteration. On option, the variable 
can be redefined to 3 6 indicating 
that simulated data will be stored 
and used. 

All logicals are initialized to FALSE in the BLOCK DATA subprogram except STAPSW. 
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/N0N2/ 



COMMON 

/NON2/XYZTP, RVTP, IOBS, INPAR, NPARAM 

Program 

Where 

Variable 

Type 

Description 

Defined 

XYZTP 

I 

A variable which defines a device 
to output XYZ type data. The 
variable is initialized to 12 for 
XYZ tape output but the ORBIT 
option card redefines the de- 
vice as the printer (6). 

BLOCK 

DATA 

OPTCRD 

IOBS 

I 

A variable which defines the 
input device for either real 
observational type data or 
rectangular coordinate type 
data. The variable is initialized 
to 10. 

BLOCK 

DATA 

OPTCRD 

INPAR 

I 

A variable reflecting whether the 
state, drag or solar radiation 
parameters are being estimated 

BLOCK 

DATA 

OPTCRD 

NPARAM 

I 

The total number of parameters 
to be estimated 

OPTCRD 


Program 

Where 

Used 


MAIN 

STATRD 


MAIN 


MAIN 

ESTIM 

READGP 

WTBMAT 

ORBIT 

PREDCT 

MAIN 

ESTIM 

SOLVGP 

WTBMAT 
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/NON2/ (continued) 


NOTE 

INPAR is initialized to 6 which indicates that state 
parameters will be estimated. The COEFGP option 
card redefines INPAR to be 0 if state parameters are 
not being estimated. Therefore, 


Label 

Value If Present 

Description 

INPAR = COEFGP 

6 

Estimate x 0 , y 0 , z 0) x 0 , y 0 , z 0 

+ADDR 

1 

Estimate C D 

+SRAD 

1 

Estimate C R 

NPARAM = INPAR 

6, 7, or 8 


+NPRTL 

NP 

Number of GP to be estimated 

+(3* NSTEST) 

NS 

Number of station coordinates 
to be estimated 

+NBIAS 

NB 

Number of biases to be 
estimated 


where 


= state parameters 


C D = drag coefficient 

C R = solar radiation coefficient 

CP = geopotential coefficients. 
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/NON3/ 


COMMON /NON3/TIMING, TTL, CDNAME, ATYPE, UNITS, BLANK, NCARDS 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

TIMING 

D 

The title "TIMING" 

BLOCK 

data 

MAIN 

STATRD 

TTL (50) 

D 

Title array for parameters 
being estimated 

BLOCK 

DATA 

MAIN 

STATRD 

CDNAME (3 2) 

D 

Array of option card names. 
There are 32 option cards 
available. 

BLOCK 

DATA 

OPTCRD 

STATRD 

ATYPE (21) 

D 

Array of data type names. 
There are 21 data types 
available. 

BLOCK 

DATA 

MAIN 

OPTCRD 

STATRD 

UNITS(15) 

D 

An array of the names of the 
units used for print-out 

BLOCK 

DATA 

MAIN 

OPTCRD 

STATRD 

BLANK 

D 

1 blank space 

BLOCK 

DATA 

MAIN 

OPTCRD 

STATRD 

NCARDS 

I 

Total number of option cards 
available (presently 32 option 

BLOCK 

DATA 

OPTCRD 


cards are available) 
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/NON4/ 


COMMON /NON4/BSEND, BSTRT, RMSTOT, EDITN, BYT PE, BSTANO, 
IRSUPR, NBIAS, NSTA, ISTEST, NSTEST, NOPRPR, IDSAT 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

BSEND(44) 

D 

End time of biases being 
estimated 

OPTCRD 

MAIN 

STATRD 

BSTRT(44) 

D 

Start times of biases being 
estimated 

OPTCRD 

MAIN 

STATRD 

RMSTOT 

R 

Input value for the desired total 
RMS. The nominal value is 
200.0. 

BLOCK 

DATA 

OPTCRD 

MAIN 

EDITN 

R 

Input sigma multiplier. The 
nominal value is 3.5. 

BLOCK 

DATA 

OPTCRD 

MAIN 

BYTPE (44) 

I 

Bias types to be estimated 

OPTCRD 

MAIN 

BSTANO (44) 

I 

Station numbers for which 
biases are requested 

OPTCRD 

MAIN 

IRSUPR(4) 

I 

Indicates which iterations to 
suppress residual printout. 
This is initialized to 0. 

BLOCK 

DATA 

OPTCRD 

STATRD 

NBIAS 

I 

Number of biases being esti- 
mated. This variable is 
initialized to 0. 

BLOCK 

DATA 

OPTCRD 

MAIN 

NSTA 

I 

Number of stations to be used 
in a run (not allowed to ex- 
ceed 50). This variable is 
initialized to 0. 

BLOCK 

DATA 

OPTCRD 

MAIN 

STATRD 

ISTEST 

I 

Array of station numbers cor- 

OPTCRD 

STATRD 


responding to those stations 
that are to be estimated 
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/N0N4/ (continued) 


Variable 

Type 

Description 

Program 

Where 

Defined 

NSTEST 

I 

Number of stations that are being 
estimated. This variable is 
initialized to 0. 

BLOCK 

DATA 

OPTCRD 

NOPRPR 

I 

Number of stations requesting 
data to be preprocessed. This 
variable is initialized to 0. 

BLOCK 

DATA 

OPTCRD 

IDSAT 

I 

Satellite identification number 

OPTCRD 


Program 

Where 

Used 


MAIN 

STATRD 

PREDCT 

ESTIM 

STATRD 


MAIN 
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/NON5/ 


COMMON /NON5/DORBIT, DAYEND, DORB1, DORB1E, EPSEC, DRATE, 
ORBRT, RATE, SEC, IEPHM, IHM, IYMD, IYBEG, IEPYMB, IYREF 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

V/here 

Used 

DORBIT 

D 

Integration start time which is 
initialized to 0.0 

BLOCK 

DATA 

OPTCRD 

MAIN 

STATRD 

DAYEND 

D 

Integration end time which is 
initialized to 0.0 

BLOCK 

DATA 

OPTCRD 

MAIN 

STATRD 

DORB1 

D 

Integration start time when 
creating an ORB1 tape is 
initialized to -1.0 

BLOCK 

DATA 

OPTCRD 

MAIN 

STATRD 

DORBlE 

D 

Integration end time for creating 
an ORB1 tape 

OPTCRD 

MAIN 

STATRD 

EPSEC 

D 

Number of seconds into the 
epoch day 

OPTCRD 

MAIN 

STATRD 

DRATE 

D 

Print time interval in days for 
a ground track request. The 
value is initialized to 999.0. 

BLOCK 

DATA 

OPTCRD 

MAIN 

STATRD 

ORBRT 

D 

Output time interval in seconds 
for an ORB1 tape 

OPTCRD 

MAIN 

STATRD 

RATE 

R 

Print time interval in seconds 
for an orbit generator run 

OPTCRD 

MAIN 

STATRD 

SEC 

R 

Number of seconds into the 
epoch day 

OPTCRD 

MAIN 

STATRD 

IEPHM 

I 

Epoch date (packed hours 
and minutes) 

OPTCRD 

MAIN 

STATRD 

IHM 

I 

Selected end time (packed 
hours and minutes) of arc 
for data reduction run 

OPTCRD 

MAIN 

STATRD 
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/N0N5/ (continued) 





Program 

Program 




Where 

Where 

Variable 

Type 

Description 

Defined 

Used 

IYMD 

I 

Selected end time (packed year, 

OPTCRD 

MAIN 



month, and day) of arc for data 


STATRD 

IYBEG 

I 

reduction run 

OPTCRB 

MAIN 



Reference year of the epoch 


STATRD 

IEPYMD 

I 

elements 

OPTCRD 

MAIN 



Date of the epoch elements 


STATRD 

IYREF 

I 

(packed year, month, and day) 

OPTCRD 

MAIN 



Reference date of the epoch 


STATRD 


elements (packed year, month, 
and day) 
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/NON6/ 


COMMON /NON6/ORBEL, SATNME, XYZNOM, PLHNOM, HN, 
SLATN, LATDN, LATMN, LONDN, LONMN, SLONN 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

ORBEL(6) 

D 

Epoch orbital elements 

MAIN 

STATRD 

SATNME 

D 

Satellite name for printout 
purposes 

MAIN 

STATRD 

XYZNOM(10, 6) 

R 

Nominal variance -covariance of 
a station position in earth fixed 
rectangular coordinate system 

STATRD 

MAIN 

PLHNOM(10, 6) 

R 

Variance-covariance of a station 
position in the geodetic coordinate 
system 

STATRD 

MAIN 

HN(10) 

R 

Nominal station height in meters 

STATRD 

MAIN 

SLATN(IO) 

R 

Number of seconds in the nominal 
station geodetic latitude 

STATRD 

MAIN 

LATDN(IO) 

I 

Number of degrees in the nominal 
geodetic latitude 

STATRD 

MAIN 

LATMN(IO) 

I 

Number of minutes in the nominal 
station geodetic latitude 

STATRD 

MAIN 

LONDN (10) 

I 

Number of degrees in the nominal 
station east longitude 

STATRD 

MAIN 

LONMN(IO) 

I 

Number of minutes in the nominal 
station east longitude 

STATRD 

MAIN 

SLONN (10) 

R 

Number of seconds in the nominal 
station east longitude 

STATRD 

MAIN 
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/OUTFOR/ 

COMMON /OUTFOR/TITLE, CENVRG, NITER 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

TITLE(36) 

D 

Title identification information 
for a run 

MAIN 

OUTPUT 

CENVRG 

R 

Convergence criteria on RMS 

MAIN 

OUTPUT 

NITER 

I 

Upper bound on total number 
of iterations for a run 

MAIN 

OUTPUT 
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/PCES/ 


COMMON /PCES/PCESW 


Variable 


PCESW 


Type 


Description 


Program 

Where 

Defined 


Program 

Where 

Used 


L 


The switch is initialized to FALSE BLOCK 
in BLOCK DATA. The PCE option DATA 
card is used to set the switch to OPTCRD 
TRUE, which indicates that rec- 
tr angular coordinate type simulated 
data will be processed. 


MAIN 

STATRD 

PREDCT 

OUTPUT 


NOTE 

Rectangular coordinate type simulated data (PCE Data) 
is processed in a fashion simular to the ordinary data types. 
The differences occur in the input format and measurement 
partial computations, where changes were made so that rec- 
tangular coordinate "measurement" types could be processed, 
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/PCOFlT/ 


COMMON /PCOFIT/NO 





Program 

Where 

Program 

Where 

Variable 

Type 

Description 

Defined 

Used 

NO 

I 

Number of observations needed 
to calculate the coefficients of 
the polynomial used to determine 
the ephemeris quantities. The 
nominal value set in the program 
is 11. 

BLOCK 

DATA 

ephqan 
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/PREBLK/ 


COMMON /PREBLK/OBS1, OBS2, IPREPR(4, 50), RFINDX(2, 50), 
INDPRE(2, 50), VHFCHN 

NOTE 

Subprogram GEOSRD is the program utilized by 
MAIN to read records of observational data, select 
and delete records on request, and preprocess optical 
and time observational data as requested by various 
indicator variables. Part of the final observational 
data is returned through COMMON/PREBLK/. The 
remainder of the observational data is returned 
through COMMON /CGEOS/. Each record of observa- 
tional data may contain one or two observations. 


Variables 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

OBS1 

D 

1st observation contained in 
record 

GEOSRD 

PROCES 

MAIN 

PROCES 

PREDCT 

OBS2 

D 

2nd observation if the record 
contains two observations 
(Refer to indicator NMEAS in 
COMMON /GEOS/.) 

GEOSRD 

PROCES 

main 

PROCES 

PREDCT 


NOTE 

The indicator variables MTYPE and NMEAS contained 
in COMMON/CGEOS/ must be used in order to in- 
terpret OBS1 and OBS2. These are used as follows : 


MTYPE = 11 OBS1 = 
NMEAS = 2 J OBS2 = 

MTYPE = 21 OBS1 = 
NMEAS = 1 J 

MTYPE = 31 OBS1 = 
NMEAS = lj 


right ascension - radians 
declination - radians 

range - meters 
range rate - meters/sec. 
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/PREBLK/ (continued) 


NOTE (continued) 



MTYPE 

NMEAS 

!i II 

OBS1 = range rate - meters/sec 

► 

(derived from Doppler) 



MTYPE 

NMEAS 

= 5) 

= 2 J 

OBS1 = £ direction cosine - 

► 

OBS2 = m direction cosine • 

dimensionless 
- dimensionless 



MTYPE 

NMEAS 

I! II 
tSD 05 

1 OBS1 = X angle - radians 
1 OBS2 = Y angle - radians 




MTYPE 

NMEAS 

= 7 1 
= 2 j 

1 OBS1 = azimuth - radians 
f OBS2 = elevation - radians 



Variables 

Type 


Description 

Program 

Where 

Defined 

Program 

Where 

Used 

IPREPR(1, N) 

I 

Indicator for preprocessing 
optical data including provisions 
for active or passive time 
correction 

MAIN 

PROCES 

GEOSRD 

(2, N) 

I 

Indication for applying new index 
of refraction 

MAIN 

PROCES 

(3,N) 

I 

Indicator for applying constant 
timing corrections 

MAIN 

PROCES 

GEOSRD 

(4, N) 

I 

Not used 


PROCES 

GEOSRD 

RFINDX(1, N) 


New values of index of refraction 
in ppm deviation from 1, from Nth 
preprocessing option card 

MAIN 

PROCES 

GEOSRD 

INDPRE(1, N) 

I 

Station number from Nth pre- 
processing option card 

MAIN 

GEOSRD 

INDPRE(2, N) 

I 

Measurement type from Nth pre- 
processing option 

main 

GEOSRD 


The observational preprocessing information contained 
in the arrays IPREPR(4, 50), RFINDX(2, 50) and 
INDPRE(2, 50) are read from the option card PREPRO. 
This information is then interpreted by subprograms 
GEOSRD and PROCES for appropriate preprocessing. 
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/PRIORI/ 


COMMON / PRIORI/BE GYMD, ENDYMD, SFLUX(2499), MGFLUX(2499), DUMMY (10000) 

VARCOV(6, 6), STASIG(3, 3, 10), PLHSIG(3, 3, 10), ELEMIN(6), XSTA(10) YSTA(IO), 
ZSTA(IO), DRAGO, EMISS0, DRAGSG, EMISSG, BBIAS(44), BIAS0(44), BIASSG(44) 

R(10, 6, 20) is equivalenced in RK to BEGYMD, thus writing over the flux tables after initial use, 

and SUM1(50, 50) is equivalenced in MAIN to BEGYMD, thus writing over the initial integration 

starting tables after initial use. 

Program Program 

Where Where 

Variable Type Description Defined Used 

BEGYMD I Beginning date of the tables of BLOCK MAIN 

solar flux and geomagnetic DATA 

activity data in the form 
YYMMDD where 

YY = last 2 digits of year 

MM = month 

DD = day 

ENDYMD I End date of the tables of solar BLOCK MAIN 

flux and geomagnetic activity DATA 

data in the form YYMMDD 

SFLUX(2499) R Daily values of solar flux at BLOCK MAIN 

2800 Me (10.7 cm) in units of DATA 

watts/ meters 2 / cycle/ second/ 
bandwidth x 10~ 22 measured at 
Ottawa ARO, adjusted to 1 
astronomical unit from BEGYMD 
to ENDYMD 

MGFLUX(2499) R Daily values of geomagnetic BLOCK MAIN 

activity indices (dimensionless) DATA 

from BEGYMD to ENDYMD 

DUMMY* (1000) R Scratch storage array 

*DUMMY is used to increase :he total length of the PRIORI COMMON block in order to store intermediate values needed in the 
computation of the starting table (see subprogram RK). 
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/PRIORI/ (continued) 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

VARCOV(6, 6) 

R 

Not used, replaced by 
/ CVRCV/VRCOV 



STASIG(3, 3, 10) 

R 

Variance- covariance matrix 
of the a priori estimate of the 

MAIN 

MAIN 

ESTIM 


geocentric earth-fixed rec- 
tangular coordinates of the 
tracking stations to be adjusted- 
meters. (Maximum of 10 sta- 
tions allowed to be adjusted 
simultaneously. ) 


PLHSIG(3, 3, 10) R 

Variance -covariance matrix of 
the a priori estimate of the geo- 
detic latitude, longitude and height 
above ellipsoid of the tracking 
stations to be adjusted— angles in 
radians, height in meters 

MAIN 
SQUANT 
(as output 
from 

subroutine 

PLHOUT) 

MAIN 

ESTIM 

ELEMIN(6) 

i D 

A priori estimate of the geocentric 
inertial rectangular coordinates of 
the position and velocity vectors of 
the satellite at the initial epoch to 
be adjusted- meters and meters/sec 



XSTA(10) " 

YSTA(10) 

ZSTA(10) 

> D 

A priori estimate of the geocentric 
. earth-fixed rectangular coordinates 
of the tracking stations to be 
adjusted— meters 

MAIN 

ESTIM 

DRAGO 

R 

C °o 



EMISSO 

R 

C *0 



DRAGSG 

R 

a ( C °o) 



EMISSG 

R 

-K) 



BBIAS(44) 

R 

b 

MAIN 

MAIN 

ESTIM 
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/PRIORI/ (continued) 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

BIAS0(44) 

R 

b o 

MAIN 

MAIN ESTIM 

BIASSG(44) 

R 

°‘< b o) 

MAIN 

MAIN ESTIM 


where : 





CL = a priori value of the satellite’s 
drag coefficient 

C Rq = a priori value of the satellite's 
reflectivity coefficient 

o-(C Do )= standard deviation of 

erf C ) = standard deviation of C„ 

\ R o / R o 

b 0 = array containing the a priori values 
of the measurement biases or station 
timing biases to be adjusted 

b= array containing the current adjusted 
values of the measurement biases or 
station timing biases 

o-(b 0 ) = array containing the standard deviations 
of b 0 . 

Data for SFLUX and MGFLUX are defined in BLOCK DATA for the periods 

BEGYMD = 650101 (January 1, 1965) 
to 

END YMD = 680430 (April 30, 1968). 

During the initialization process of the main program, values for 15 days of solar flux and geomag- 
netic data, with the first value corresponding to the day of the epoch of the initial elements, are 
placed in another common block (COMMON/FLXBLK/). Appropriate provisions are made for esti- 
mating values if the epoch of the initial elements lies outside the range of the table. After these 
data have been placed in COMMON/FLXBLK/, the storage area in COMMON/PRIORI/, starting with 
the variable BEGYMD and ending with the last value of MGFLUX, is utilized by subroutines RK and 
ESTIM to store temporary integration arrays and the matrix SUM1(50, 50). As a result the stored 
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/PRIORI/ (continued) 


values of solar flux and geomagnetic activity, data are overlaid and no longer available for use. 
This fact prevents "stacking" cases "back to back" since the required solar flux and geomagnetic 
data are available only on the 1st pass through the program. 

After solution, the matrix SUM1(50, 50) contains the variance-covariance matrix of all parameters 
which are adjusted by the DC. 
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/RKST/ 


COMMON /RKST/RSTEP 





Program 

Where 

Program 

Where 

Variable 

Type 

Description 

Defined 

Used 

RSTEP 

D 

Starting table stepsize. It should 
generally be a fraction of the 
Cowell stepsize. The nominal 
value is 20 seconds but can be 
changed using the RKSTEP option 
card. 

BLOCK 

DATA 

OPTCRD 

RK 
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/RKT/ 


COMMON /RKT/CSTEPT 





Program 

Program 




Where 

Where 

Variable 

Type 

Description 

Defined 

Used 

CSTEPT 

D 

Required Cowell stepsize time 

TABLE 

RK 


points for the starting table 
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/SCRTCH/ 


Program 

MAIN 


SQUANT 


COEFF 


COMMON /SCRTCH/DUMMY(500) 

This common block is used as a scratch storage area by 
the following programs: 

MAIN 

COEFF 

SQUANT. 

It is only used for the storage of information that need 
not be saved. Usage of this common block in other 
routines is permitted under the same restrictions. 

Usage of /SCRTCH/ 

To store station positions in output format at the beginning of a case and at 
the end of each iteration when station position adjustment has been requested. 
Used to pass station variance-covariance information to subroutine SQUANT 
when station position adjustment has been requested.. 

To obtain station positions at the beginning of each case for conversion to 
earth-fixed rectangular coordinates. To return station positions to the 
main program, in output format, when station position adjustment has been 
requested. 

To store matrices needed for calculation of coefficients for polynomial fit. 
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/SETSW/ 


Variable 

IND(9) 


110 


INDX 


COMMON / SETSW/lND(9), 110, INDX 


Type 


Description 


Program 

Where 

Defined 


Program 

Where 

Used 


I 


I 


I 


An array of indicators which define SWTEST 
computed GO TO statements in the 
force computations for the varia- 
tional equations 

Ari indicator which when set to: SWTEST 

1- indicates state parameters will 
not be estimated, 

2- indicates that only state param- 
eters will be estimated, 

3- indicates that state and force 
model parameters are to be 
estimated. 

An indicator which when set to 0 SWTEST 
is the same as 110 set to 1; a value 
of 6 is the same as 110 set to 2 or 3. 


VEVAL 


VEVAL 

MAIN 

ESTIM 

DRAG 

FRCS 


VEVAL 

MAIN 

ESTIM 

DRAG 

FRCS 


NOTE 

Indicators 110 and INDX are used to pack the GRPAR, 
PXPXO and PMPXO arrays so that if a particular 
parameter type is not to be estimated, the lower 
parameter partials move up. 
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/SIGMAC/ 


COMMON / SIGMAC/ SIGCHG, IMTYPE, ISTNO, SIGSTD, CULL, NS1G, NCULL 


This common block is used to transfer optional input card information to the observation 
processing modules. 


Variable 

Type 

Description 



The arrays, SIGCHG, IMTYPE, ISTNO, store the optional inputs 
used to override the standard observation sigmas defined by the 
observation types in SIGSTD. This permits specification of an 
observation sigma by station and/or observation type. 

BIGCHG(50) 

R 

SIGCHG(J) = Value of sigma in units of observation type (see 
SIGSTD) 

IMTYPE (50) 

I 

Observation type for SIGCHG(J) 

ISTNO (50) 

I 

Applicable station number for SIGCHG(J) 

SIGSTD(14) 

R 

Pre-stored array of observation sigmas by type, which is 
initialized in BLOCK DATA to the following: 


Type 

Sigma Value 

Observation 

(1) 

20. seconds of arc 

rt. ascension 

(2) 

25. meters 

range 

(3) 

10. meter s/sec 

range rate 

(4) 

1 . 

Doppler 

(5) 

0.3 

direction cosine, 

(6) 

50. degrees 

X angle 

(7) 

50. seconds of arc 

azimuth 

(8) 

20. seconds of arc 

declination 

(12) 

0.3 dimensionless 

direction cosine, m 

(13) 

50. degrees 

Y angle 

(14) 

50. seconds of arc 

elevation 

(16) 

100. meters 

geocentric coordinate x 

(17) 

100. meters 

geocentric coordinate y 

(18) 

100. meters 

geocentric coordinate z 

(19) 

0.1 meters/sec 

geocentric coordinate x 

(20) 

0.1 meters/sec 

geocentric coordinate y 

(21) 

0. 1 meters/sec 

geocentric coordinate z 
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/SIGMAC/ (cdntinued) 


Vai’iable Type 


Description 


CULL(2, 100) 


NSIG 

NCULL 


I Defines the range of observations, by sequence number to be 

deleted from a DC run; where CULL(1, K) is the first number 
and CULL(2, K) is the last number. K is the number of sets of 
observations to be deleted. 

I Number of variables in arrays, SIGCHG, IMTYPE, ISTNO which 

is initialized to 0 in BLOCK DATA. 

I Number of sets of observations to be deleted, K. Initialized to 

0 in BLOCK DATA. 


These variables are defined in subprogram OPTCRD from optional input cards and are subsequently 
used by subprograms DODSRD and GEOSRD. 
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/STANUM/ 


COMMON /STANUM/NAME (50), ISTANO(50) 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

NAME(N) 

D 

Identifying name of the Nth track- 
ing station (assumed to have a 
maximum of six alphanumeric 
characters) 

BLOCK 

DATA 

MAIN 

STATRD 

PLHOUT 

ISTANO(N) 

I 

Identifying number of the Nth track- 
ing station (currently assumed to 
have a maximum of four digits) 

BLOCK 

DATA 

STATRD 

MAIN 


The arrays NAME and ISTANO are also stored by 
STATRD as the station position cards are read in. 
A maximum of 50 tracking stations are allowed by 



GEOSTAR-I. 
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/VRBLOK/ 


COMMON /VRBLOK/RTXYSQ, COSLAM(31), SINLAM(31), RR 


Variable 

Type 

Description 


Program 

Where 

Defined 

Program 

Where 

Used 

RTXYSQ 

D 

y'x 2 + y 2 - meters 


EGRAV 

EGRAV 

veval 

COSLAM (N) 

R 

cos (N - l)x; N = 1, 2, • 

• • 31 

EGRAV 

BLOCK 

DATA 

EGRAy 

VEVAL 

SINLAM(N) 

R 

sin (N - 1)\; N = 1, 2, • 

•• 31 

EGRAV 

BLOCK 

DATA 

EGRAV 

VEVAL 

RR 

D 

r - meters 


EGRAV 
(same 
value as 
R placed 
in 

COMMON 

/XYZ/) 

VEVAL 


where: 


x, y = geocentric inertial x and y 
components of the position 
vector of the satellite 

r = geocentric distance of 
satellite 

k = geodetic longitude ox sub- 
satellite point. 


BLOCK DATA initially defines: RTXYSQ =0 

COSLAM(l) = 1 
COSa-AM( 2) - COSLAM (31) = 0 
SINLAM(l) - SINLAM(31) = 0 
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/VRBLOK/ (continued) 


Computing Note 


Particular attention should be paid to the fact that 
subroutine VEVAL establishes the variables in this 
block as: 

COMMON /VRBLOK/A1, G2, CSLM(31), SNLM(30), R1 
DOUBLE PRECISION Al, Rl 

The correspondence of the /VRBLOK/ variables in 
VEVAL with the /VRBLOK./ variables in EGRAV are 

EGRAV 


RTXYSQ 

COSLAM(l) 

COSLAM(2) - COSLAM (31) 
SINLAM(l) 

SINLAM(2) - SINLAM(3 1) 
RR 


■VEVAL 

Al 

G2 

CSLM(l) - CSLM(30) 
CSLM(3 1) 

SNLM(l) - SNLM(30) 
Rl 



/WORKER/ 


COMMON /WORKER/X(20, 3, 50), XD(20, 3, 50), XDD(20, 3, 50), 
XXDD(3, 50), CDEL(2), TT(2), T(2), KI(2), N(2), NEQ, IPA(2) 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

X(20, 3, 50) 

D 

Array of position vectors for the 

ORBIT 

CSTEP 

XD(20, 3, 50) 

D 

equations of motion and the 
variational equations in meters 

Corresponding array of velocity 

ORBIT 

SUMS 

HEMINT 

RK 

FRCS 

TABLE 

TABLE B 

CKDIFF 

CSTEP 

XDD(20, 3, 50) 

D 

vectors in meters/second 
Corresponding array of accelera- 

ORBIT 

SUMS 

HEMINT 

RK 

FRCS 

TABLE 

TABLE B 

CKDIFF 

CSTEP 

XXDD(3, 50) 

L 

tion vectors in meters/second 2 
Current acceleration vector at a 

ORBIT 

SUMS 

HEMINT 

RK 

FRCS 

TABLE 

TABLES 

CKDIFF 

CSTEP 



given time point for the equations 
of motion and the variational 
equations 


SUMS 

HEMINT 

RK 

FRCS 

TABLE 

TABLEB 

CKDIFF 
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/WORKER/ (continued) 





Program 

Program 




Where 

Where 

Variable 

Type 

Description 

Defined 

Used 

CDEL(2) 

D 

CDEL(l) = h x 

MAIN 

CSTEP 



CDEL(2) = h 2 

OPTCRD 

SUMS 




BLOCK 

HEMINT 




DATA 

RK 





FRCS 





TABLE 





TABLEB 





CKDIFF 

TT(2) 

D 

TT(1) = h/ 

ORBIT 

CSTEP 



TT(2) = h 2 2 


SUMS 





HEMINT 





RK 





FRCS 





TABLE 





TABLEB 





CKDIFF 

T(2) 

D 

T(l) is the time of the current 

ORBIT 

CSTEP 



position, velocity, and acceleration 


SUMS 



vectors from epoch for the equa- 


HEMINT 



tions of motion. 


RK 



T(2) is the time of the current 


FRCS 



position, velocity, and acceleration 


TABLE 



vectors from epoch for the varia- 


TABLEB 



tional equations. 


CKDIFF 

KI(2) 

I 

KI(1) indicates the position (in the 

ORBIT 

CSTEP 



position, velocity, and acceleration 


SUMS 



arrays) of the current position, 


HEMINT 



velocity and acceleration vectors 


RK 



for the equations of motion. 


FRCS 



KI(2) indicates the position (in the 


TABLE 



position, velocity, and acceleration 


TABLEB 



arrays) of the current position, 


CIIDIFF 



velocity, and acceleration vectors 





for the variational equations 
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/WORKER/ (continued) 





Program 

Pr )gram 




Where 

Where 

Variable 

Type 

Description 

Defined 

Used 

N(2) 

I 

N(l) = P l - 2 

ORBIT 

CSTEP 



N(2) = P 2 - 2 


SUMS 

HEMINT 

RK 

FRCS 

TABLE 

TABLEB 

CKDIFF 

NEQ 

I 

Number of equations being 

MAIN 

CSTEP 



integrated 


SUMS 

HEMINT 

RK 

FRCS 

TABLE 

TABLEB 

CKDIFF 

IPA(2) 

I 

IPA(l) is the storage location of 

ORBIT 

CSTEP 



the interpolated vectors for the 


SUMS 



equations of motion; 


HEMINT 



^A(2) is the storage location of 


RK 



the interpolated vectors for the 


FRCS 



variational equations 


TABLE 

TABLEB 

CKDIFF 


where 

hj = stepsize used to integrate the equations of motion 
h 2 = stepsize used to integrate the variational equations 
p l = order of formulas used to integrate the equations of motion 
p 2 = order of formulas used to integrate the variational equations. 
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/XYZ/ 


COMMON /XYZ/X, Y, Z, XDOT, YDOT, ZDOT, R, RSQ, RQ, TI 





Program 

Program 




Where 

Where 

Variable Type 

Description 

Defined 

Used 

xl 

1 

Satellite position vector in 

FRCS 

EGRAV 

Y 

r D 

meters 

(Defined 

DRAG 

Z J 

i 


from vari- 
ables input 
through 
COMMON 
/WORKER/ 
from CSTEF 
or RK 

VEVAL 

xdotI 


Satellite velocity vector in 

FRCS 

EGRAV 

YDOT 

► D 

meters/second. 

(Defined 

DRAG 

zdotJ 



from vari- 
ables input 
through 
COMMON 
/'WORKER/ 
from CSTEP 
or RK 

VEVAL 

R 

D 

r - meters 

EGRAV 

DENSTY 




FRCS 

DRAG 




VEVAL 

EGRAV 

VEVAL 

ESQ 

D 

r 2 - meters 2 

FRCS 

VEVAL 

VEVAL 

RQ 

D 

r 3 - meters 3 

FRCS 

VEVAL 

VEVAL 

TI 

D 

Current integration time of 

FRCS 

VEVAL 



the variational equations. 

CSTEP 
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/XYZ/ (continued) 


where : 


!x\ 

Y 

w 

/x\ 

Y 

W 


geocentric inertial rectangular coordinates of the po- 
sition and velocity vectors of the satellite at the epoch 
when subroutine FRCS is called by either subroutine 
CSTEP or RK for the calculation of accelerations 
during the integration process, r and v are passed by 
CSTEP and RK to FRCS by means of the COMMON 
/WORKER/. If the equations of motion and the vari- 
ational equations are being integrated with the same 
stepsize, then FRCS places this information in COMMON 
/XYZ/ for use in subprograms DRAG, EGRAV, DENSTY, 
and VEVAL. However, if two different stepsizes are be- 
ing used, then f and v are interpolated vectors at the time 
needed by the variational equations and placed in the 
COMMON /XYZ/ by CSTEP for use in VEVAL. 
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/XYZOUT/ 


COMMON /XYZOUT/XYZEND(6) 

NOTE 

The array XYZEND(6) contains the geocentric inertial 
rectangular coordinates of the position and velocity 
vectors of the satellite at each observation time or 
requested print time. 





Program 

Program 




Where 

Where 

Variable 

Type 

Description 

Defined 

Used 

XYZEND(l) 

D 

x - meters 

MAIN 

MAIN 




(as out- 

(as in- 

XYZEND(2) 

D 

y - meters 

put from 

put to 

XYZEND(3 ) 

D 

z - meters 

subroutine 

subroutine 




ORBIT) 

ORBIT) 





PREDCT 

XYZEND(4) 

D 

x - meters/second 



XYZEND(5) 

D 

y - meters/second 



XYZEND(6) 

D 

z - meters/second 




where 





r = 


v = 



geocentric inertial rectangular coordinates of the 
position and velocity vectors of the satellite 
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5.2 NONAME -GEOSTAR-I ODP COMMON Block Cross Reference Table 

This section contains a cross reference table describing the COMMON area structure in the 
GEOSTAR-I ODP. 
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COMMON SLOCKS 


GEOSTAR-I COMMON Cross Reference Table 


SUBROUTINES 
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GEQSTAR-I COMMON Block Cross' Reference Table (Continued) 
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Nofe... Common Blocks COWS, ERR22, and SSTEP are used by the Variable Step Program 
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5.3 LUNGFISH— GEOSTAR-I SOLVE COMMON Blocks 

This section contains a detailed description of the COMMON areas used in the GEOSTAR-I 
SOLVE program. 
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/EARTH/ 


COMMON /EARTH/AE, FLAT 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

AE 

D 

Semi-major axis of earth 

MAIN 

OPSTAT 

FLAT 

D 

Flattening coefficient 

MAIN 

OPSTAT 
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/INPUT/ 


COMMON /INPUT/IOPT1, IOPT2, IOPT3, MEORG, MEPAR, MESUP, 
MERED, MEDIC, MECOM, MEINV, NCMSUP, JCMSUP(500), 
IDIN, BNAMIN(3), NBSUP, JBSUP(500), NARG 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

IOPT1 

I 

This variable can have the values: 

(0) — check input data cards and 

B matrix tape; 

(1) — check input cards only; 

(2) — no checks will be made, 

for illegal values or format 

MAIN 

INCK 

IOPT2 

I 

This variable can have the values : 

(0) — normal solve run 

(1) — process only one input 

B matrix 

MAIN 

INCH 

IOPT3 

I 

This variable can have the values: 

(0) — intermediary matrices gen- 

erated during a solution 
saved on tape #30. 

(1) — intermediary matrix to be 

printed 

MAIN 

INCK 

MEORG 

I 

Indicates the type of edit requested 
for the original input B matrix 

MAIN 

INCK 

MEPAR 

I 

Indicates the type of edit requested 
for the original input parameter 
set matrices 

MAIN 

INCK 

MESUP 

I 

Indicates the type of edit requested 
for the suppressed matrices 

MAIN 

INCK 

MERED 

I 

Indicates the type of edit requested 
for the reduced matrices 

MAIN 

INCK 

MEDIC 

I 

Indicates the type of edit requested 
for the backsubstitution matrices 

MAIN 

INCK 
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/INPUT/ (continued) 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

MECOM 

I 

Indicates the type of edit requested 
for the final combined matrix 

MAIN 

INCK 

MEINV 

I 

Indicates the type of edit requested 
for the inverse of the final com- 
bined matrix 

MAIN 

INCK 

NCMSUP 

I 

The number of parameters to be 
suppressed from the combined 
matrix 

MAIN 

MAIN 

INCK 

JCMSUP (500) 

I 

Array of parameter labels to be 
suppressed from the combined 
matrix 

MAIN 

MAIN 

INCK 

IDIN 

I 

Identification number of each 
B matrix used 

MAIN 

MAIN 

INCK 

BNAMEN(3) 

R 

B matrix name 

MAIN 

MAIN 

NBSUP 

I 

Number of parameters to be sup- 
pressed from each B matrix 

MAIN 

MAIN 

INCK 

JBSUP(500) 

I 

Array of parameter labels to be 
suppressed from the current 
B matrix 

MAIN 

MAIN 

INCK 

NARG 

I 

Total number of parameters to 
be suppressed from the current 
B matrix 

MAIN 

MAIN 


NOTE 

Before calling SUPRSS and LBLSUP, the total list of 
parameters to be suppressed is formed from the 
JBSUP and JCMSUP arrays; likewise for the variable 
NARG. 
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/LOCAL/ 


LOCAL is a common area used in each subroutine to reduce storage requirements of the pro- 
gram. Data transfer occurs only between the routines INVERT, MINV, EDIT, OP ARC, OPGRAV, 
OPSTAT, where the array A(200, 200), a double precision array containing the matrix to be in- 
verted and the inverted matrix, is transferred. 
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/PERM/ 


COMMON /PERM/B(500), BNEW(500), SIGTON, VI, V2, V3, VIC, V2C, V3C, IRT1, IDMAT, 
NROW, NCOL, NOB, ITYPE, BNAME(3), IRT2, LABS(500), LABSI(500), 
NRONI, LABSN(500), NROWN, NGRAV, NSTAT, NARC, JPAGE, LINE, 
NOREC, NOREC2, DBVS, ALPHA(8), IDCOMB, NR^WC, NCOLC, NOBC, 
ITYPE C, BNAMEC(3), LABSC(501) 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

B(500) 

D 

The right hand side of the 
normal equations 

INVERT 

MINV 

BNEW(500) 

D 

The parameter set after 
backsubstitution 

BACKSB 

BACKSB 

SIGTON 

D 

The predicted variance 

MAIN 

MAIN 

VI, V2, V3 1 

VIC, V2C, V3C J 

D 

Arc variances and 
total variances after com- 
bining arcs 

MAIN 

UPCOMB 

INVERT 

ELIM 

MAIN 

MAIN 

IRTl 

I 

Integer indicating the first 
record in a given matrix 



IDMAT 

I 

The ID number of a given 
matrix; on the header record 
of a matrix tape 



NROW 

I 

The number of rows in a given 
matrix; on the header record 
of a matrix tape 



NCOL 

I 

The number of columns in a given 
matrix; on the header record of 
a matrix tape 



NOB 

I 

The number of observations for 
each arc, and also for the com- 
bined arcs 
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/PERM/ (continued) 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

ITYPE 

I 

Matrix type identifier and used 
to indicate edit type 



BNAME(3) 

R 

An array containing the name of 
the B matrix being processed 



IRT2 

D 

Integer indicating the second 
record in a given matrix 



LABS (500) 

I 

An array containing the 
parameter labels 

MAIN 

UPCOMB 

SUPRSS 

LBLSUP 

LABSI(500) 

I 

An array containing the 
parameter labels of the 
inverted matrix 

INVERT 

BACKSB 

NROWI 

I 

Number of rows in the 
inverted matrix 

INVERT 

BACKSB 

LABSN(500) 

I 

An array containing the 
parameter labels of the 
backsubstitution matrix 

BACKSB 

MAIN 

NROWN 

I 

Number of rows in the back- 
substitution matrix 

BACKS B 

MAIN 

NGRAV 

I 

Number of geopotential 
coefficient parameters 

CALTYP 


NSTAT 

I 

Number of station position 
parameters 

CALTYP 


NARC 

I 

Number of arc dependent 
parameters 

CALTYP 


JPAGE 

I 

Page number of printed output 

CHECK 


LINE 

I 

Line number of printed output 

CHECK 
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/PERM/ (continued) 





Program 

Program 




Where 

Where 

Variable 

Type 

Description 

Defined 

Used 

NOREC, NOREC2 

I 

Integers used to control back- 
spacing of tapes 

MAIN 

BEDIT 

OBVS 

R 

Set equal to NOB to compute 
sigma 

MAIN 

MAIN 

ALPHA(8) 

R 

An array containing alpha- 
numeric run identification 
printed at the head of each 
page of output 

MAIN 

CHECK 

IDCOMB 

I 

ID number of the combined 

MAIN 

MAIN 



matrix 


COMB 

INCK 

NROWC 

I 

Number of rows in combined 
matrix 

UPCOMB 

COMB 

NCOLC 

I 

Number of columns in com- 
bined matrix 

UPCGMB 

COMB 

NOBC 

I 

The total number of observa- 
tions over all arcs used 

UPCOMB 

COMB 

ITYPEC 

I 

The ITYPE value for the out- 
put combined matrix 

MAIN 

COMB 

BNAMEC(3) 

R 

Output name of the combined 
matrix 

MAIN 

COMB 

LABSC(501) 

I 

The label array for the com- 

UPCOMB 

ELIM 



bined matrix 


COMB 
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/SOGMA/ 


COMMON /SOGMA/XA(500) 





Program 

Where 

Program 

Where 

Variable 


Description 

Defined 

Used 

XA(500) 

D 

An array containing the stand- 
ard deviations of the parameters 

INVERT 

OPARC 

OPSTAT 

OPGRAV 
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5.4 LUNGFISH— GEOSTAR-I SOLVE COMMON Block Cross Reference Table 

The following table details the COMMON area structure of the GEOSTAR-I SOLVE program. 


SUBROUTINES 
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VI. PROGRAM FLOW CHARTS AND SUBROUTINE SUMMARIES 


The following sections contain the flow diagrams of the MAIN or control routines in the 
GEOSTAR-I ODP and SOLVE programs. A summary of all the modules used in these programs 
and cross reference tables detailing the subroutine structure is also presented. 
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6.1 GEOSTAR-I ODP Flow Diagrams 

This section contains detailed flow diagrams of the GEOSTAR-I ODP executive program 
MAIN and the orbit generator control subprogram ORBIT. 
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6.2 GEOSTAR-I ODP Subroutine Summary 

This section contains a brief summary of all modules used in the GEOSTAR-I ODP. Modules 
of the original NONAME system are included as well as those developed for GEOSTAR-I ODP. 
Complete documentation of the NONAME modules can be found in Reference 1. 
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ADDYMD 

AND2 

APPER 

BLKSTA 

CLEAR 

CKDIFF* 

COEFF 

COEFL 

CSTEP 

DARCTN 

DATES 

DAYEAR 

DENMUL 


Adds or subtracts an integral number of days from a date packed into 
one word in the form YYMMDD giving a new packed date. 

A function which determines the logical "AND" of two arguments, each 
two bytes long. The FORTRAN function LAND is called. 

Computes apogee and perigee heights in kilometers given the Keplerian 
elements, radius of the earth, and inverse of flattening. 

Sets up the working arrays of station positions merging initial values, 
pre stored via the BLOCK DATA subprogram with card input station 
positions. The arrays are ordered with stations to be estimated first . 

Stores zeros in any two dimensional REAL*4 array. 

Computes the kth backward difference which is needed to estimate 
the local truncation error in Runge-Kutta generated points. 

Estimates coefficients of a polynomial of a requested degree (maxi- 
mum 7th degree) by the method of least squares using a maximum of 
20 observations. 

Prints non- zero coefficients of the spherical harmonic expansion of 
the geopotential model used by the system. 

Solves a set of 1st and 2nd order, ordinary differential equations 
using the summed Stormer/Adams-Bashforth predictor and the 
Cowell/Adams -Moulton corrector. 

Function returns the arc-tangent of (Y/X) in radians, between 0 
and 2v. The FORTRAN function DATAN2 is called. 

Converts days elapsed from Jan. 0.0 of an input reference year into 
a three word date of the form YYMMDD, HHMM, SEC. 

Converts days elapsed from Jan. 0.0 of the reference year into two 
words of integral number of days and integral number of seconds. 

Evaluates a 3rd degree polynomial given the coefficients of the poly- 
nomial and a value of the independent variable. 


•Subroutines used only in the variable stepsize version of GEOSTAR-1 ODP. 
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DENSTY 

DIFF 

DINRAD 

DJUL 

DNVERT 

DODSRD 
DOT PRD 
DRAG . 
EGRAV 

ELEM 

ELEMK 

EPHQAN 


A function which computes atmospheric density dependent on height 
and temperature. Temperature is derived from the Jacchia-Nicolet 
model considering: solar activity, semiannual variation, diurnal bulge, 
and geomagnetic activity. 

Calculates difference between any two dates in the 20th centry, (input- 
dates input in two words of the form YYMMDD, HHMMSS; output- 
integral days and seconds of a day). 

Converts angles in arc measurements or time measurements to 
radians (input: integer degrees or hours, integer minutes, REAL *8 
seconds). 

Converts input days since Jan. 0.0 of reference year to an adjusted 
Julian date (reference Jan. 0, 1950). 

Double precision matrix inversion using Gauss-Jordan method of 
condensation with partial (column) pivoting. No restrictions on di- 
mension of matrix. 

Reads and processes observation data in the DODS system format to 
set up the system working arrays of input observations. 

A function which computes the dot product of two three-dimensional 
vectors. 

Computes accelerations in rectangular coordinates on a satellite due 
to drag forces. 

Computes acceleration in rectangular coordinates on a satellite due 
to geopotential forces, and acceleration partials (forcing functions) 
in rectangular coordinates due to geopotential coefficients. 

Converts inertial position and velocity vector to osculating orbital 
elements. Inputs and outputs are in COMMON. 

Converts orbital elements as ELEM with input and output in the 
argument list. 

Executive routine calling various subprograms which evaluate the 
three components of the moon's inertial unit vector and geocentric 
distance, the sun's inertial unit vector and geocentric distance, and 
the equation of equinoxes. 
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EQN 

EQXJATR 

ERROR 

ESTIM 

FRCS 

FIT 

GDET 

GEOSRD 

GTIMIN, GTIMOT 
HEMINT 

INV2/INV3 


Computes nutation in longitude, obliquity and right ascension; true 
obliquity of date. 

Transforms the rectangular coordinates of a vector referenced to 
either the mean or true equator and equinox of one epoch to either 
the mean or true equator and equinox of another epoch. The subroutines, 
PRECES and NUTATE are called. 

Provides printed messages for error conditions. 

Sums all measurement partials into the normal equations matrix and 
right hand sides for each observation data point. When all observations 
have been processed, a priori parameter weights (input to the system 
as parameter sigmas or a full weight matrix) are applied. 

Executive routine calling various subprograms which evaluate ac- 
celerations in rectangular coordinates on a satellite due to the vari- 
ous forces. 

Evaluates polynomial functions given the coefficient matrix and 
a value of the independent variable. There is no restriction 
on the maximum degree or number of the polynomials to be 
evaluated. 

Evaluates the determinant of a matrix (reduction to diagonal form 
using elementary row and column operations). There is no restric- 
tion on the dimension of the matrix. 

Reads and processes input observation data in the NASA Science Data 
Center format and sets up the working arrays of observation data 
points. Data may be selected by measurement type, and/or station 
and/or time period. 

Measures the elapsed time of each iteration and computes the total 
elapsed time of a run in hundredths of seconds. 

Interpolates for position, velocity and position partials at a time 
between equally spaced points, using the Hermite interpolation 
formula. 

Inserts a matrix by using the Gauss-Jordan method of condensation 
with partial (column) pivoting. 


189 



MMATRX Multiplies an n X n matrix with an n x m matrix where n need not 

equal m. 

MULMAT Multiplies two input matrices. 


NUMBER 


NUTATE 


OBSDOT 


Searches the entries of members of an array and compares with an 
input number or bit configuration. Index number of the array member 
whose entry is matched is returned. Zero is returned if no match is 
found. It is employed as a device to avoid excessive use of DO loops. 

Generates the nutation angles to transform from mean equator and 
equinox to true equator and equinox. 

Calculates time derivatives of requested observation types. It is 
used in the computation of time biases. 


OPT CRD 


ORB1 


ORBIT 


OUTPUT 

OUTRAD 

PLHOUT 

POSVEL 

PRECES 


Reads cards which redefine stored values of earth, satellite, and 
integration parameters as well as initializing various program options. 

Provides control to generate an ephemeris tape in the ORB1 format. 
ORBIT is called for the ephemeris points. 

Executive program W'hich receives a state vector and its epoch, 
initializes required constants, and utilizes an integrator subprogram 
and an interpolation subprogram to find the state vector at a new 
epoch, as well as the position partials of any physical parameters 
being estimated. 

Prints values of the earth, satellite, and integration parameters. 

Converts radians to degrees or hours, minutes and seconds. 

Converts a tracking station location and variance- covariance matrix 
in geodetic rectangular coordinates to geodetic latitude, longitude and 
height coordinates. 

Converts osculating orbital elements to inertial position and velocity 
vectors. Input and output is in COMMON. 
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Generates the angles for precession from mean equator and equinox 
of one epoch to mean equator and equinox of another. The year 1950 
is used as a base year. 



PREDCT Converts inertial position and velocity to a computed observation re- 

quired for the residuals (O-C). The partials of an observation with 
respect to state or parameters at epoch are computed as: 

JM_ _ <9M dx t 

3x 0 dx t 5x 0 

where: 

M is an input observation measurement at time t, 

x represents state or parameters to be estimated in the DC, 

where the subscript zero denotes epoch time, 

j~~ is computed in PREDCT, 

dx t 

is computed in ORBIT. 

PRNTPR Provides printed information when preprocessing of input observa- 

tion data is requested. 

PROCESS Preprocesses a requested observation type. 

READGP Reads the input data cards dnd sets up working arrays for the geo- 

potential estimation problem, 

REFCOR Transforms vectors to mean equator and equinox of reference year. 

The subroutine PRECES is called. 

REFIMP Computes nominal ionospheric refraction correction for range and 

range rate data. 

RK Solves a set of 1st order ordinary differential equations using an 8th 

order Runge-Kutta integration method. 

ROTMAT Generates the rotation matrix given the rotation angle and axis. 

RYMDI Separates date packed into one word in the form YYMMDD into the 

three words YY, MM, and DD. 

SATCLC Utilizes known differences in time between the GEOS-I satellite clock 

and UTC, and, assuming that the input time is the recorded UTC time 
of an observed active flash observation, calculates the time difference 
between the UTC of the actual flash and the input UTC. 

SATCL2 Serves the same purpose as SATCLC for GEOS-2. 
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SDOLL Used with NASA Data Center format tapes to provide for the selec- 

tion of input observation data. Reads up to 100 SELECT or DELETE 
cards, then reads records of observational data from an input tape or 
disk and writes these records onto a second tape or disk unit, select- 
ing or deleting records by time interval, and/ or station and/ or type of 
observation as determined by the SELECT or DELETE cards. 

SDTPRD Computes the single precision dot product of two three-dimensional 

vectors. 


SOLVGP 


SQUANT 


STAINF 

STATRD 

STORGP 

SUMS 

SUMTOB 


Computes the parameter corrections using matrix inversion given 
the weighted normal equations as produced by ESTIM. A gradient 
option is provided as an alternate method if normal matrix inverson 
is inadequate. 

Used in the initialization phase of a DC run to convert geodetic spheri- 
cal coordinates (geodetic latitude, longitude and height above compu- 
tational spheroid) to geodetic rectangular coordinates of a tracking 
station, and components of the station Zenith, East and North unit 
vectors for an array of stations (maximum 50). Also computes on first 
pass matrix of partial derivatives of geodetic rectangular coordinates 
with respect to geodetic spherical coordinates for tracking stations 
whose coordinates are to be adjusted (maximum 10), On subsequent 
calls to the program, geodetic rectangular coordinates of the stations 
being adjusted together with their variance-covariance matrices are 
transformed to geodetic spherical coordinates. Their local East, 

Zenith and North unit vectors are recomputed on the basis of the 
adjusted positions. 

Computes statistical information at the end of each DC iteration; in- 
cluded are residual summaries by station and data tape. 

Redefines and prints both the geodetic spherical coordinates and the 
geodetic rectangular coordinates of the stations (maximum 50). 

Sets up the working arrays of the geopotential model used by VEVAL. 

Computes che first and second sums necessary for the summed form 
of the predictor -corrector multistep integration formulas. 

Provides for the shift of a specified matrix row and array variable 
in the subprogram WTBMAT. 
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SUN 

SUNGRV 

SWTEST 

SYMMET 

TABLE 

TABLEB* 

TDIF 

TEST* 

VCONV 


VEVAL 


Computes the solar position unit vector components and radius 
vector and the mean equinox and ecliptic of date. 

Computes the acceleration in rectangular coordinates due to solar 
and lunar gravity. 

Computes 'indices which identify the forces being applied to the equa- 
tions of motion and the forces being applied to the variational equations. 

Computes the elements of given square matrix below the main 
diagonal on the assumption that the matrix is symmetrical. 

Executive routine calling other subprograms to produce the initial 
table of starting values for the Cowell integrator. 

Produces a new table of starting values for the Cowell integrator 
whenever a stepsize change occurs. 

A function that computes differences between the time systems Al, 
UTC, UT2 and UT1. 

Computes the local truncation error and tests this error against a 
set of tolerances to determine if an increase or decrease of the in- 
tegration stepsize is necessary. If a stepsize change occurs, this 
subprogram computes the new stepsize and starting values by calling 
TABLEB. 

Converts a variance-covariance, VIN, matrix from one system to 
another, VOUT, by computing the matrix product B = P T AP from 
input matrices A and P where A and P are both 3x3. It is used 
for station information where : 

B = VOUT, in spherical coordinates 
A = VIN, in rectangular coordinates 

P = Partials of VOUT variables with respect to VIN variables. 

Computes the partials of a satellite's acceleration vector with re- 
spect to inertial position and velocity vectors and physical parameters 
to be estimated. The force model includes geopotential harmonics 
through 4th order zonal and 3rd order tesseral harmonics, drag, 
lunar gravity, solar gravity, and solar radiation. 


‘Subroutines used only in the variable stepsize version of GEOSTAR-I ODP. 
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WTBMAT 


XEFIX, YE FIX, 
XINERT, YINERT 


YMDAY 


Used in multi-arc runs to reorder the normal equations and right 
hand sides, and to provide a parameter list and labels in the B matrix 
format (see Appendix B). 

Functions that convert earth-fixed rectangular coordinates to in- 
ertial rectangular coordinates and vice versa, specifically: 

XEFIX - given inertial X and Y return earth fixed X 

YE FIX - given inertial X and Y return earth fixed Y 

XINERT - given earth fixed X and Y return inertial X 

YINERT - given earth fixed X and Y return inertial Y. 

Computes days elapsed from Jan. 0.0 of a reference year to input 
date in form YYMMDD, HHMM, SEC, Reference year is set by pro- 
gram to be equal to year of input data on first call to program. 
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6.3 GEOSTAR-I ODP Subroutine Cross Reference Tables 

This section contains cross reference tables detailing the subroutine structure in the 
GEOSTAR-I ODP executive program MAIN and the orbit generation control subprogram ORBIT. 
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CALLED SUBROUTINES 


MAIN Program Cross Reference Table 


ADDYMD 

AND2 

APPER 

COEFF 

COEFL 

DARCTN 

DATES 

DAYEAR 

DENMUL 

DENSTY 

DIFF 

DINRAD 

DJUL 

DODSRD 

QQTPRD 

DRAG 

EGRAV 

ELEM 

ELEMK 

EQN 

EQUATR 

ERROR 

ESTIM 

FRCS 

GDET 

GEOSRD 

GTIMIN 

GTIMOT 

HEMINT 

!NV2 

INV3 

MMATRX 

MOONAD 

MULMAT 

NUMBER 

NUTATE 

OBSDOT 

OPTCRD 

ORB1 

ORBIT 

OUTPUT 

POSVEL 

PRECES 

PREDCT 

READGP 

REFCOR 

ROTMAT 

RYMDI 

SATCL2 

SATCLC 

SDOLL 

SDTPRD 

SOLVGP 

SQUANT 

STAINF 

STATRD 

STORGP 

SUN 

SUNGRV 

SYMMET 

TDIF 

VEVAL 

WTBMAT 

XEFIX 

XINERT 

YEFIX 

YINERT 

YMDAY 


CALLING SUBROUTINES 

MBCCDDDDDDDEE EFGNOOO 
A LOSAAE IJORLP QRE UBPr 
I.KETTYNFUDAehUCOTSTb 

N SFEEESFLS G MQASSADC 1 

tfpsat R at R T O " 

A RY D NRDET 


ISmBSi 
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CALLED SUBROUTINES 


MAIN Program Cross Reference Table (Continued) 
CALLING SUBROUTINES 



O 

p 

p 

p 

p 

R 

R 

R 

S 

s 

s 

s 

s 

s 

T 

T 

T 

T 

V 

w 

Y 


R 

L 

R 

R 

R 

E 

E 

K 

A 

A 

O 

Q 

T 

T 

A 

A 

D 

E 

E 

T 

M 


B 

H 

E 

E 

O 

F 

F 

T 

T 

L 

U 

A 

A 

B 

B 1 

s 

V 

B 

D 


1 

O 

C 

D 

C 

C 1 

C 

C 

V 

A 1 

T 

L 

L 

F 

T 

A 

M 

A 


T 

U 

E 

C 

/ E 

O 

M 

L 

L 

G 

N 

N 

R 

E 

E 

L 

A 

Y 



T 

S 

T 

's 

R 

P 

2 

C 

P 

T 

F 

D 

B 

T 

BLKSTA 














B 








CKDIFF 


















B 




CLEAR 













B 









COEFL 














B 








CSTEP 

El 





















DARCTN 




B 


















DATES 














B 








DENSTY 



















B 



DIFF 









B 












B 

DINRAD 









■ 



B 










DNVRT1 











B 











DOTPRD 




B 

B 














B 



EGRAV 



















B 



EPHQAN 

D 





















FIT 




Q 


















FRCS 

B 







B 








a 






HEMINT 

B 















B 






MULMAT 



B 



B 
















NUMBER 






1 

B 







B 








NUTATE 






B 
















OBSDOT 





B 














m 

B 


OUTRAD 












B 







m 

m 


PLHOUT 












B 


B 








PRECES 






B 






B 










PRNTPR 














B 








PROCES 


■ 


B 


















REFCOR 

B 


















B 



REFIMP 





B 

















RK 















B 







ROTMAT 



B 



















SQUANT 














B 

n 







SUMS 

B 















B 






SUMTOB 




















B 


SWTEST 

B 

■ 




















SYMMET 











a 









B 


TABLE 

B 





















TABLEB 














mi 

m 



B 




TEST 

B 













n 

5 







VCONV 


m 












B 








VEVAL 

B 







B 






LJ 


X 






XSF1X 




X 











m 







XINERT 




X 










■ 

■ 

n 






YEFIX 




X 


















YINERT 




X 


















YMDAY 



a 

- 





D 

B 







D 






Note... Subroutines CKDIFF, TABLEB, and TEST are used by the Variable Step Program 
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CALLED SUBROUTINES 


ORBIT Subprogram Cross Reference Table 


CALLING SUBROUTINES 

OCCDDDDE FNPR RTTTVY 

ROSEI JRPRUREKAAEEM 

betnfuahctef bbsvd 

ifesflgqsacc lltaa 

tfpt a teo ee ly 

Y N E s R B 


CKDIFF 

COEFF 

CSTEP 

DENMUL 

DENSTY 

DIFF 

DJUL 

DOTPRD 

DRAG 

EGRAV 

EPHQAN 

EQN 

FRCS 

HEMINT 

INV2 

INV3 

MMATRX 

MOONAD 

MULMAT 

NUTATE 

PRECES 

REFCOR 

RK 

ROTMAT 

RYMDI 

SUMS 

SUN 

SUNGRV 

SWTEST 

TABLE 

TABLEB 

TDIF 

TEST 

VEVAL 

YMDAY 



Note... Subroutine CKDIFF, TABLEB, and TEST are used by the Variable Step Program 
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6.4 GEOSTAR-I SOLVE Flow Diagram 

This section contains detailed flow diagrams of the GEOSTAR-I SOLVE executive programs 
in the normal execution mode and in the "invert only" mode. 
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SOLVE (NORMAL RUN) 


THIS TAPE CAN BE RE- 
INTRODUCED AS INPUT IN 
A LATER RUN IF A RELATED 
CB MATRIX AND PARAMETER 
SET ARE INTRODUCED AS 
WELL. 

USE TAPE 29 


THIS TAPE CAN BE RE- 
INTRODUCED AS INPUT IN 
A LATER RUN JF A RELATED 
CB MATRIX AND D''C 
MATRIX ARE INTRODUCED 
AS WELL. 

USE TAPE 19 



( RETURN ) 


♦NOTE: THE COMBINED ARCS RUN DOES NOT ELIMINATE ANY STATE PARAMETERS, NOR BACKSUBSTITUTE TO RETRIEVE 

THEM. IT IS THUS IN ITS LOGIC FLOW CLOSER TO THE "NORMAL" SOLVE RUN, WHILE CLOSER IN ITS ACTIVITIES TO THE 
INVERT ONLY RUN. THE COMBINED ARCS RUN BORROWS FROM BOTH. 

THE FINAL COMBINED ARCS MATRIX HAS NO CORRESPONDING D"’c MATRICES, AND HAS ITS PARAMETER SET WRITTEN 
OUT AFTER THE MATRIX ON THE SAME TAPE. HENCE, TO REINTRODUCE THE COMBINED ARCS MATRIX INTO A LATER RUN, 
TREAT IT AS A SELF-CONTAINED B MATRIX. 
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UP- 

PRESSED \ I ETER 

MATRIX, j \ MATRIX 

TAPE 28 J \ TAPE 19 











6.5 GEOSTAR-I SOLVE Subroutine Summary 


This section contains a brief summary of all the modules used in the GEOSTAR-I SOLVE. 
Modules of the original LUNGFISH system are included as well as those developed for the 
GEOSTAR-I system. A more complete documentation of the LUNGFISH modules can be found in 
Reference 2. 
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ANDREE 

Computes the Penrose pseudoinverse of a given matrix. The com- 
putational rank of the given matrix is also determined. 

BACKSB 

Solves the backsubstitution equations for the arc dependent param- 
eters using the backsubstitution matrices generated during the 
elimination process, the associated right hand sides, and the arc 
independent parameter solution set. 

BEDIT 

Writes the matrices' to be edited onto unit 30 together with an edit 
format code. 

CALTYP 

Calculates the number of gravity coefficient, station position and 
arc dependent parameter labels in a given SOLVE run. 

CHECK 

Checks the number of lines printed out on any given page, and if 
enough lines have been written, a new page is begun, headed by the 
alphanumeric input of data card 1. 

COMB 

Performs the matrix combination of the reduced matrices resulting 
from arc dependent elimination, and writes the resultant combined 
matrix on unit 28. 

DARCTN 

Function returns the arc tangent of Y/X in radians between 0 and 
27 t . The FORTRAN function DATAN2 is called. 

EDIT 

Prints the matrices to be edited (stored on unit 30) in the required 
edit format. 

ELIM 

Performs a Gauss-Jordan elimination of the arc dependent parameters 
from a given matrix, producing a reduced matrix containing only arc 
independent parameters and a backsubstitution matrix for the eliminated 
parameters. 

ERROR 

Prints an error message identifying a particular error condition and 
terminates the run. 

GNORM 

Normalizes a given gravity coefficient. 

INCK 

Checks the format of the various input options, values, and 
matrices. If these formats are not correct, ERROR is called. 

INVERT 

Control routine for matrix inversion. Calls for either inversion or 
pseudoinversion of a given matrix. 

203 



LBLSUP 

MATSUP 

MINV 

OPARC 

OPSTAT 

OUTRAD 

PHLINN 

SORTX 

SUPRSS 


Sets the labels corresponding to suppressed parameters to zero. 

Performs the matrix suppression by deleting rows and 'iumns and 
writes the resulting suppressed matrix on unit "X 

Inverts and solves a matrix equation using a Gauss-Jordan elimination 
method, with a pivot search on the diagonal elements. 

Prints and punches the new estimates of the arc dependent parameters. 



Prints and punches the new estimates of the station position param- 
eters. Output is in both geodetic rectangular and spherical coordinates. 

Converts radians to degrees, minutes and seconds. 

Converts station coordinates from geodetic rectangular coordinates 
to geodetic latitude, longitude and height above spheroid. 

Sorts the gravity parameter labels and values so that all the C nm 
geopotential parameters are in ascending order immediately preceding 
the S nm geopotential parameters. 

Control program for matrix suppression. 


UP COMB 


Produces the combined matrix identification and label records. 
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6.6 GEOSTAR-I SOLVE Subroutine Cross Reference Tables 

The following table details the subroutine structure of the GEOSTAR-I SOLVE program. 
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CALLED SUBROUTINES 


CALLING SUBROUTINES 



M 

B 

B 

c 

E 

E 1 

1 

M 

o 

o 

o 

p 

s 

u 


A 

A 

E 

o 

D 

L 

N 

N 1 

p 

p 

p 

H 

u 

p 


1 

C 

D 

M 1 1 

C 

V 
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A 

s 

G 
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p 

c 


N 

K 1 a 

T 

M 

K 

E 

V 

R 

T 

R 1 

R 

o 



s 

T 

R 

C 

A 

A 

N 

s 

M 



B 

T 


T 

N 

V 

s 

B 















" 


ANDREE 


Hi 






X 









BACKSB 

X 

HI 











mm 



BEDIT 

mm 

mm 






mm 

Hi 




m 



CALTYP 

mm 

mm 


M 



X 

m 

I 

mm 

X 

X 


X 

X 

CHECK 

mm 




X 


X 


X 

X 

X 

mm 

H 



COMB 

mm 




■ ■ 








B 



DARCTN 





■ 








X 



EDIT 

mm 



H 












ELIM 

mm 















ERROR 

mm 

X 

mm 

X 

X 

X 

X 



mm 






GNORM 



B§gf 









X 




INCK 

X 


I 





■ 








INVERT 

X 









mm 






LBLSUP 










B 




X 


MATSU P 



■ ■ 








H 



mm 


MINV 
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X 



X 



B 





OPARC 

X 


B 







B 






OPGRAV 

mm 















OPSTAT 

X 















OUTRAD 


■ 









mm 





PLHINN 


■ ■ 


■ ■ 

s 



■ ■ 



X 





SORTX 




■ ■ 

Hi 



■ 




X 




SUPRSS 
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UPCOMB 
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Appendix A 

Variable Stepsize Version of GEOSTAR-I ODP 
New and Updated Modules in GEOSTAR-I ODP 

This section describes the m edifications and additions made to the subroutine structures of the 
GEOSTAR-I ODP for the variable stepsize version of the GEOSTAR-I ODP program. 

To obtain the variable stepsize version, the following GEOSTAR-I ODP modules were modified: 


MAIN 

ORBIT 

FRCS 

OPT CRD 

CSTEP 

HEMINT 

OUTPUT 

TABLE 

RK 

BLOCK DATA 


SUMS 

ERROR 


SWTEST 


In addition, the following modules were developed specifically for the variable stepsize 
GEOSTAR-I ODP: 

CKDIFF 
TABLE B 
TEST 

It should be noted that this modified version of GEOSTAR-I ODP has all the capabilities of the 
original GEOSTAR-I ODP, including the multiple arc capability. The single exception is that the 
modified program can only solve for a total of 20 unknowns, as opposed to 50 for the original 
version. 

In the following sections, those modules listed above -which are new, or significantly modified 
existing subroutines will be documented in detail. The remaining modules which received only 
minor modifications are briefly outlined, with changes indicated. 
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A.l Modifications to Existing GEOSTAR-I ODP Modules 


The modifications made to the GEOSTAR-I ODP are designed to- 


• Call the new subroutines (see Sec. 6.1 Orbit Flow) 

• Extend computational algorithm to allow for stepsize modification in the integration of the 
equations of motion and the variational equations 

• Extend tables to provide space for storing position', velocity and position partial vectors to 
allow for stepsize modification. 


A summary of these modules, as modified for the variable stepsize GEOSTAR-I ODP, follows: 


• MAIN 

FRCS 

HEMINT 

• RX 
SUMS 
SWTEST ^ 

OPT CRD 1 

• OUTPUT > 
ERROR J 

• INPUT 


• ORBIT 


• TABLE 


• CSTEP 


The ODP control program. Modified to determine when to use the 
variable stepsize option. 


Modified to extend the position, velocity, and position partial vectors 
tables to allow for stepsize modification. 


Modified to include the I/O required for the stepsize modification 
process. 

Modified to extend the available program options to include variable 
stepsize and new data base constants required for this new option. 

The integration control program. Modified to control the new sub- 
routine TEST as well as the variable stepsize option logic. 

The control program for the initial table of starting values to be used 
in the Cowell integrator. Modified to halve or double the stepsize of 
both the equations of motion and the variational equations if neces- 
sary. In addition, print statements have been included to indicate 
a stepsize change. 

The subroutine that solves a set of 1st and 2nd order differential equa- 
tions using the predictor -corrector technique. Modified to save the 
predicted and corrected values of position which are used for local 
error computation. 
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A. 2 New Variable Stepsize ODP Modules 
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CKDIFF 


Purpose; 

To compute the kth backward difference of the accelerations which is needed to estimate the 
local truncation error of the initial starting table for the equations of motion. 

Called By: 

TEST 

Method: 

To compute the kth difference of a table of function values, the formula 

k . \ 

V*f(x) = 22 f ( x " ih ) 

i= 0 


is used. 

Calling Sequence: 


CALL CKDIFF (DIFF) 


COMMON Blocks Used: 



WORKER 


Variables Not in COMMON: 

FORTRAN Name 

Format 

Description 

DIFF (3) 

D 

Kth backward difference vector 

BINC(20) 

D 

Array of binomial coefficients 
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TEST 


Purpose : 

To compute and test the local truncation error to determine whether an increase or a decrease 
of the stepsize s of the equations of motion and the variational equations is required to satisfy the 
local error constraint equation. 

Called By: 

ORBIT 

TABLE 

Calls : 

TABLEB 

CKDIFF 

Method: 


The subroutine computes the local error Rn and determines whether the constraint equation 

T0L2 < Rn < TOL1 

is satisfied, where TOL1 and TOL2 are specified error bounds. If this equation is satisfied, control 
returns to the calling program. If Rn > TOLl, the stepsize is decreased and if Rn > TOL2, the step- 
size is increased. The methods used to compute the local error Rn, the new stepsize and the new 
tabular points depend on the calling program. Let k be the number of backpoints to be used; then 
we have: 

Case I: Calling Program TABLE 

In this case the test is made to determine whether the initial Cowell stepsize used in the 
starting table is adequate. The local error is estimated by the formula 

Kn = V k X k 


where the kth difference V k is computed by calling subroutine CKDIFF. The new stepsize com- 
puted is either half or double the original stepsize, and the new table of points are computed by 
returning to TABLE and restarting with the newly computed stepsize. 
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Case II: Calling Program ORBIT 


In this case the test is made to determine whether the current Cowell stepsize at the nth in- 
tegration step is adequate. The local error is estimated by the formula 

Kn ar C|x n p -x n c | , 

where c is an input error constant and x n p , x n c are the predicted and corrected values of the satellite 
position vectors. The new stepsize is computed by the formula 

/TOL3\ 1/k 

^ncw ^old \ Rn / 


where TOL3 is a specified value for the allowable local error, TOL2 < TOL3 < TOL1. In this case 
the new table of points at this new stepsize is computed by calling TABLEB and interpolating for 
the necessary values at the new stepsize. 

The nominal values for the tolerances used in the program are: 

TOL1 = .25 x 10" 4 meters 
TOL2 = .25 x 10" 10 meters 
TOL3 = ,25 x 10" 7 meters 


which can be changed on option. 
Calling Sequence 


CALL TEST 


COMMON Blocks Used: 


ERR22 WORKER 

COWS LIMITS 

SSTEP 


Variables Not in COMMON: 
FORTRAN Name 


Format 

D 


Description 
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DIFF(3) 

TP 


D 


Kth backward difference 

Time span available at d 
stepsize 



FORTRAN Name 

Format 

Description 

TPP 

D 

Time span needed at new 



stepsize 

OSTEP 

D 

Computed new stepsize 

FN1 

D 

K 


Error Constants 

The error constants for the equations of motion when the calling program is ORBIT are: 


Order 

Value 

4 

0.00000000 

5 

.050000000 

6 

.052631579 

7 

.048721340 

8 

.044032445 

9 

.039713131 

10 

.035961476 

11 

.032748940 

12 

.029998465 

13 

.027632051 

14 

.025582503 

15 

.023794802 


The error constants for the equations of motion when the calling program is TABLE are: 

Value 

.0041666666 
.0041666666 
.0036541005 
.0031415344 
.0027086089 
.0023553241 
.0020677822 
.0018320857 
.0016369383 
.0014736450 
.0013356018 
.0012177854 


Order 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 
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TABLEB 


Purpose: 

To produce a new starting table for both the equations of motion and the variational equations 
whenever a stepsize change occurs. 

Called By: 

TEST 


Calls: 

HEMINT 

FRCS 

SUMS 

VEVAL 

Method: 


The subroutine produces a new set of position, velocity, and position partials for the equa- 
tions of motion and the variational equations respectively at the new stepsize by: 

(i) Computing the times at which points are to be produced; 

(ii) Calling HEMINT with these times to interpolate for the necessary starting values. 

A new starting table of acceleration vectors for the equations of motion and the variational 
equations is then obtained by using the newly computed position and velocity vectors in the FRCS 
and VEVAL subroutines. Next, new first and second sums for both the equations of motion and the 
variational equations using position, velocity, partials, and acceleration vectors at the new step- 
size are computed by calling SUMS. 

This subroutine also prints the value of the new stepsize and the time from epoch when the 
stepsize change occurred. 


Calling Sequence : 


COMMON Blocks Used: 


CALL TABLEB 


WORKER GRBLOK 
LIMITS ANPART 
SSTEP 
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Variables Not in COMMON: 


FORTRAN Name 

Format 

Description 

TPP 

D 

Time span needed at new 
stepsize 

TMP 

D 

Time span for interpola- 
tion of points 

TIME 

D 

Time of last good point 
computed 

TI 

D 

Interpolation time 

M 

I 

Interpolation order 

KIP 

I 

Lower index of points to 
be interpolated 

K1 

I 

Starting index of pants to 
be interpolated 

XS(15, 3, 20) ^ 
XDS(15, 3, 20) j 

f 

Storage arrays for inter- 
polated position, velocity 
and position partials at 
the new stepsize 


217 



A. 3 Variable Stepsize COMMON Block Variable Description 

This section contains a detailed description of the COMMON areas used in the GEOSTAR-I 
ODP variable o .psize program which supplement those given in Section 5.1. 
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/cows/ 


COMMON /COWS/TOL1, TOL2, TOL3, STPMIN, MODE 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

TOL1 

D 

Upper bound on local error 

BLOCK 
DATA 
OPT CRD 

TEST 

TOL2 

D 

Lower bound on local error 

BLOCK 
DATA 
OPT CRD 

TEST 

TOL3 

D 

Desired local error at new 
integration stepsize 

BLOCK 

DATA 

OPTCRD 

TEST 

STPMIN 

D 

Minimum stepsize allowed 

BLOCK 

DATA 

OPTCRD 

TEST 

MODE 

I 

Variable indicating the mode 
of operation of the program. 
The values of this variable 
are: 

MAIN 

BLOCK 

DATA 

OPTCRD 

ORBIT 

TABLE 


(2) Variable step mode 

(4) Fixed step mode. 

The nominal values of these variables built into the program are: 

TOL1 = ,25 x 10" 2 * 4 * * meters 

TOL2 = .25 x 10" 10 meters 

TOL3 - .25 x 10" 7 * * meters 

STPMIN = 5.0 sec 

MODE = 2. 


219 



/ERR22/ 


COMMON /ERR22/SPO(3), SPOO(3) 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

SPO(3) 

D 

Saved predicted value of posi- 
tion vector for local error test 

CSTEP 

TEST 

SPOO(3) 

D 

Saved corrected value of posi- 
tion vector for local error test 

CSTEP 

TEST 
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/SSTEP/ 


COMMON /SSTEP/TEMP3, NNS, NORDER 


Variable 

Type 

Description 

Program 

Where 

Defined 

Program 

Where 

Used 

TEMP3 

D 

Old stepsize when a, stepsize 
change occurs 

TEST 

TABLE B 

NNS 

I 

Saved N(l) = Pj - 2 

TEST 

TABLE B 

NORDER 

I 

Saved ORDER(l) 

TEST 

TABLEB 


where : 

Pj = order of the formula used to integrate the equations of motion 
ORDER (1) = Pj 
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Appendix B 

B Matrix Tape Format 


Record 

No. 

1 

2 

3 

4 


N + 2 
N + 3 
N +4 
N + 5 
N + 6 


Record 

Size 

60 

4N + 8 
8N + 12 
8N + 12 
8N + 12 

4N +4 
8N +4 
60 


Description of Contents 

Header — matrix size (N), total variance (VI), etc. 
Matrix labels 

Matrix data - row 1 (starting with right hand side) 

Matrix data - row 2 (starting with right hand side) 

Matrix data - row N (starting with right hand side) 

Matrix parameter set identification 

Matrix parameter labels 

Matrix parameter values 

End of logical tape 


RECEDING PAGE BLANK NOT FILMrr, 
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Record No. 1 — Header 


Item 

Format 

No. of Bytes 

Value/Limits 

Record Type 

1*4 

4 

10001 

Matrix Identification Number 

I * 4 

4 

1 to 99998 

No. of Matrix Rows, N 
(= No. of Parameters) 

I * 4 

4 

N < 200 

No. of Matrix Elements Per Row 
(Including the Right Hand Side) 

l * 4 

4 

1 + N 5 201 

Total Variance, VI 

R * 8 

8 


Weighted Variance, V2 

CO 

* 

pa 

8 


Arc Variance, V3 

R * 8 

8 


No. of Observations 

I * 4 

4 


Matrix Type 

1*4 

4 

8 

Matrix Name (alphanumeric) 

R * 4 

12 

12 EBCDIC 
characters 


Total No, of Bytes 

60 
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Record No. 2 — Matrix Labels 


Item 

Format 

No. of Bytes 

Value/Limits 

Record Type 

I * 4 

4 

10002 

Dummy 

1*4 

4 

0 

Parameter Labels 

1*4 

4N 



Total No. of Bytes 

4N + 8 



Record No. (J + 2) — Matrix Data — 

Row J (J = 1, 2, 

- • • , N) 

Item Format 

No. of Bytes 

Value/ Limits 

Record Type 1*4 

4 

10003 

Elements of Jth Matrix Row R * 8 

(Starting with the Right Hand 

Side) 

8(1 + N) 

N 200 




Total No. of Bytes 8N + 12 
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Record No. (N f 3) — Matrix Parameter Set Identification 


Item 

Format 

No. of Bytes 

Value/Liraits 

Record Type 

1*4 

4 

10011 

Matrix Identification No.* 

1*4 

4 

1 to 99998 

Dummy 

1*4 

4 

1 

No. of Parameters* 

1*4 

4 

N < 200 

Dummy 

I * 4 

28 

All 0 

Code for Parameter Set 

I * 4 

4 

8 

Matrix Name* 

R * 4 

12 

12 alphameric 
characters 


Total No. of Bytes 60 


*Same as in Record No. 1 . 
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Record No. (N + 4) — Matrix Parameter Labels 


Item 

Format 

No. of Bytes 

Value/ Limits 

Record Type 

1*4 

4 

10012 

Parameter Labels* 

1*4 

4N 

N < 200 


Total No. of Bytes 4N + 4 



Record No. (N + 5) — Matrix Parameter Values 

Item Format No. of Bytes 

Record Type 1*4 4 

Parameter Values R * 8 8N 

Total No. of Bytes 8N + 4 


Value/Limits 
10013 
N i 200 


Item 
Record Type 
Dummy 


Record No. (N + 8) — End of Logical Tape 


Format 


I * 4 
I * 4 

Total No. of Bytes 


No. of Bytes 
4 
56 


60 


Value/Limits 
-19991 
All 0 


227 



